*" σ面上で PV (ちゃんと 3 次元)を求める 96/09/02 akahori ********************************************************************* PROGRAM GTPV2 IMPLICIT NONE * #ifdef SYS_IBMS INCLUDE (GTSINC) INCLUDE (GZSIZE) #else #include "gtsinc.F" #include "gzsize.F" #endif COMMON /GMWORK/ MWORK REAL MWORK ( IJKDIM ) * *"[入出力データの領域確保] CHARACTER HTHETA ( NDC )*(NCC) !" θ REAL GTHETA ( IJKDIM ) CHARACTER HU ( NDC )*(NCC) !" ｕ REAL GU ( IJKDIM ) CHARACTER HV ( NDC )*(NCC) !" ｖ REAL GV ( IJKDIM ) CHARACTER HPS ( NDC )*(NCC) !" Ps REAL GPS ( IJKDIM ) CHARACTER HPV ( NDC )*(NCC) !" Ps REAL GPV ( IJKDIM ) * *"[途中の計算に必要なワーク領域] CHARACTER HUS ( NDC )*(NCC) !" ∂ｕ∂σ REAL GUS ( IJKDIM ) CHARACTER HVS ( NDC )*(NCC) !" ∂ｖ/∂σ REAL GVS ( IJKDIM ) CHARACTER HVOR ( NDC )*(NCC) !" ∂ζ/∂σ REAL GVOR ( IJKDIM ) CHARACTER HTHETAX ( NDC )*(NCC) !" ∂θ/∂ｘ REAL GTHETAX ( IJKDIM ) CHARACTER HTHETAY ( NDC )*(NCC) !" ∂θ/∂ｙ REAL GTHETAY ( IJKDIM ) CHARACTER HTHETAS ( NDC )*(NCC) !" ∂θ/∂σ REAL GTHETAS ( IJKDIM ) * *"[データに初期値を与える] INTEGER ITHETA, IU, IV, IPS, JPV DATA ITHETA / 50 / DATA IU / 52 / DATA IV / 54 / DATA IPS / 56 / DATA JPV / 64 / * CHARACTER THT *(NFILN) DATA THT / '$GTTMPDIR/THETA' / CHARACTER U *(NFILN) DATA U / '$GTTMPDIR/U' / CHARACTER V *(NFILN) DATA V / '$GTTMPDIR/V' / CHARACTER PS *(NFILN) DATA PS / '$GTTMPDIR/PS' / CHARACTER PV *(NFILN) DATA PV / '$GTTMPDIR/PV' / * *"[球面調和関数による微分のための定義] LOGICAL DX !" T: x 微分 (DY=Tの方が優先される) LOGICAL DY !" T: y 微分 LOGICAL A !" T: 地球半径で割る INTEGER WB !"(cosφ)**WB の重みを微分前に掛ける (x微分の場合WBと INTEGER WA !"(cosφ)**WA の重みを微分後に掛ける WAは同じ働き) LOGICAL ZONAL !" T: 一番目の軸を y 軸とみなす(波数０成分を仮定) * CHARACTER ITEM *(NCC) DATA ITEM / 'PV' / CHARACTER UNIT *(NCC) DATA UNIT / 'PVU' / CHARACTER TITLE *(NCC*2) DATA TITLE / 'PV' / CHARACTER DSET *(NCC) DATA DSET / ' ' / CHARACTER EDIT *(NCC) DATA EDIT / ' ' / CHARACTER ETTL *(NCC) DATA ETTL / ' ' / LOGICAL APND DATA APND / .FALSE. / LOGICAL GRESET DATA GRESET / .TRUE. / LOGICAL HELP DATA HELP / .FALSE. / * EXTERNAL GSPLIN2, CORIOL, DSPHE, U2VOR, PTLVOR * INTEGER NFILE, NOPT, IOS, IAX, IEOTHT, IEOU, IEOPS, NDIF, & IXDIM, IXSTR, IXEND, & IYDIM, IYSTR, IYEND, & IZDIM, IZSTR, IZEND, & IEOV NAMELIST /OPTION/ & THT, U, V, PS, !" 入力ファイル名 & PV, !" 出力ファイル名 & ITEM, UNIT, TITLE, DSET, EDIT, ETTL, & APND ,GRESET, HELP *"[オプション] CALL OPTARG ( 91, 'OPTION', 'HFILE', NOPT, NFILE ) READ (91,OPTION,IOSTAT=IOS) CLOSE(91) IF ( IOS.NE.0 .OR. HELP ) THEN WRITE(6,OPTION) STOP ENDIF *"[gtool をはじめるよ] CALL GTOPEN *"[データ入力のための前処理] CALL GTSIZE ( HTHETA , IJKDIM ) CALL GTSIZE ( HU , IJKDIM ) CALL GTSIZE ( HV , IJKDIM ) CALL GTSIZE ( HPS , IJKDIM ) CALL GTSIZE ( HPV , IJKDIM ) CALL GTSIZE ( HUS , IJKDIM ) CALL GTSIZE ( HVS , IJKDIM ) CALL GTSIZE ( HVOR , IJKDIM ) CALL GTSIZE ( HTHETAX, IJKDIM ) CALL GTSIZE ( HTHETAY, IJKDIM ) CALL GTSIZE ( HTHETAS, IJKDIM ) CALL GMSIZE ( IJKDIM ) *"[データ入力] CALL GFROPN ( ITHETA , THT ) CALL GFROPN ( IU , U ) CALL GFROPN ( IV , V ) CALL GFROPN ( IPS , PS ) CALL GFOOPN ( JPV , PV , APND ) *"[鉛直補間のためのパラメータ設定] IAX=3 !" 3 番目の軸を補間 CALL SETPAR2 ( IAX ) * 1100 CONTINUE * *"[データ入力] CALL GFREAD !" θ 入力 O ( HTHETA , GTHETA , IEOTHT , I ITHETA , 1 ) * CALL GFREAD !" Ｕ 入力 O ( HU , GU , IEOU , I IU , 1 ) * CALL GFREAD !" Ｖ 入力 O ( HV , GV , IEOV , I IV , 1 ) * CALL GFREAD !" Ps 入力 O ( HPS , GPS , IEOPS , I IPS , 1 ) * IF ( ( IEOTHT .EQ. 0 ) .AND. ( IEOU .EQ. 0 ) + .AND. ( IEOV .EQ. 0 ) .AND. ( IEOPS .EQ. 0 ) ) THEN * *"[θの鉛直微分] CALL GPFSET + ( HTHETA , GTHETA , + ' ' , ' ' , + HTHETAS, GTHETAS ) NDIF=1 !" 一階微分 CALL SETDIF ( NDIF ) CALL GMCAL1 !" dθ/dσ を求める I ( GSPLIN2, I HTHETAS, GTHETAS, I ' ' , ' ' ) *"[uの鉛直微分] CALL GPFSET + ( HU , GU , + ' ' , ' ' , + HUS , GUS ) NDIF=1 !" 一階微分 CALL SETDIF ( NDIF ) CALL GMCAL1 !" dｕ/dσ を求める I ( GSPLIN2, I HUS , GUS , I ' ' , ' ' ) *"[vの鉛直微分] CALL GPFSET + ( HV , GV , + ' ' , ' ' , + HVS , GVS ) NDIF=1 !" 一階微分 CALL SETDIF ( NDIF ) CALL GMCAL1 !" dｖ/dσ を求める I ( GSPLIN2, I HVS , GVS , I ' ' , ' ' ) *"[θの X 微分] CALL GPFSET + ( HTHETA , GTHETA , + ' ' ,' ' , + HTHETAX, GTHETAX ) DX=.TRUE. DY=.FALSE. ZONAL=.FALSE. WB=0 WA=-1 A=.TRUE. CALL EDSPHE(DX,DY,WA,WB,A,ZONAL) CALL GMCAL1 I ( DSPHE , M HTHETAX, GTHETAX, I ' ' , ' ' ) *"[θの Y 微分] CALL GPFSET + ( HTHETA , GTHETA , + ' ' ,' ' , + HTHETAY, GTHETAY ) DX=.FALSE. DY=.TRUE. ZONAL=.FALSE. WB=0 WA=0 A=.TRUE. CALL EDSPHE(DX,DY,WA,WB,A,ZONAL) CALL GMCAL1 I ( DSPHE , M HTHETAY, GTHETAY, I ' ' , ' ' ) *"[相対渦度の鉛直成分] CALL GPFSET + ( HU , GU , + ' ' ,' ' , + HVOR , GVOR ) CALL GMCAL2 I ( U2VOR, M HVOR , GVOR, !" input U ; output VOR I HV , GV , I ' ' , ' ' ) *"[相対渦度-->絶対渦度] CALL GMCAL1 !" 絶対渦度を求める I ( CORIOL, M HVOR , GVOR, I ' ' , ' ' ) *"[合体] CALL GUQSIZ I ( HVOR , O IXSTR , IXEND, IXDIM, O IYSTR , IYEND, IYDIM, O IZSTR , IZEND, IZDIM ) CALL GHCOPY ( HVOR , HPV ) CALL PTLVOR I ( GPS , I GUS , I GVS , I GVOR , I GTHETAS, I GTHETAX, I GTHETAY, O GPV , I IXDIM, IYDIM, IZDIM ) *"[ PV データ出力] IF ( ITEM .NE. ' ' ) THEN CALL GHCSET( HPV , 'ITEM', ITEM ) ENDIF IF ( UNIT .NE. ' ' ) THEN CALL GHCSET( HPV , 'UNIT', UNIT ) ENDIF IF ( TITLE .NE. ' ' ) THEN CALL GHCSTS( HPV , 'TITL', TITLE ) ENDIF IF ( DSET .NE. ' ' ) THEN CALL GHCSET( HPV , 'DSET', DSET ) ENDIF IF ( GRESET ) THEN CALL GHRSGP( HPV ) ENDIF IF ( GRESET ) CALL GHRSGP( HPV ) CALL GFWRIT I ( HPV , GPV , I JPV , 1 , 0 ) * GOTO 1100 ENDIF * commented out by toyoda * CALL GTCLSE STOP END ******************************************************************** *" (旧 95/01/16 horinout CUBIC SPLINE INTERPOLATIO) *" 改良している *" 補間か微分かを決める NDIF は ENTRY 文で代入 *=================================================================== SUBROUTINE GSPLIN2 I ( HHEAD1, GDATA1, O HHEADO, GDATAO, D IDIM1 , JDIM1 , KDIM1 ) IMPLICIT NONE * #ifdef SYS_IBMS INCLUDE (GZSIZE) #else #include "gzsize.F" #endif * INTEGER IDIM1, JDIM1, KDIM1 CHARACTER HHEAD1 ( * )*(*) REAL GDATA1 ( IDIM1,JDIM1,KDIM1 ) CHARACTER HHEADO ( * )*(*) REAL GDATAO ( * ) * internal work CHARACTER HAX1 ( NDC )*(NCC) INTEGER NAXD PARAMETER (NAXD=2048) REAL AX1(NAXD) REAL X(NAXD),Y(NAXD),Y2(NAXD) * INTEGER IAX, NDIF, IAXE, NDIFE SAVE IAX, NDIF INTEGER I, J, IJKDO, N, K, IJK, IEOD1 REAL VMISS, AX, YB * * / read axis / * CALL GTSIZE(HAX1,NAXD) CALL GUQAXV !" original axis I ( HHEAD1, IAX, 'LOC' , O HAX1, AX1 , IEOD1 ) * * / header / * CALL GHPGET( HHEAD1,'MISS',VMISS ) * IJKDO=IDIM1*JDIM1*KDIM1 CALL GTSIZE( HHEADO,IJKDO ) * * / cubic spline interpolation / * DO 1300 J=1,JDIM1 DO 1300 I=1,IDIM1 N=0 DO 1310 K=1,KDIM1 IF(GDATA1(I,J,K).NE.VMISS) THEN N=N+1 X(N)=AX1(K) Y(N)=GDATA1(I,J,K) ENDIF 1310 CONTINUE CALL SPLINE(X,Y,N,1E31,1E31,Y2) DO 1320 K=1,KDIM1 IJK=I+IDIM1*(J-1)+IDIM1*JDIM1*(K-1) AX=AX1(K) CALL SPLINT(X,Y,Y2,N,AX,NDIF,YB) GDATAO(IJK)=YB 1320 CONTINUE 1300 CONTINUE * RETURN *=========== ENTRY SETPAR2(IAXE) IAX = IAXE RETURN *=========== ENTRY SETDIF (NDIFE) NDIF=NDIFE RETURN END ************************************************************************ *(Original Source Location: /home/horinout/src/numrecipes/splin.f) * * NUMERICAL RECIPES * *" 3.3 3次スプライン補間 *" subroutines SPLINE 準備 ( Y2 = y" ) *" SPLINT 補間 仮定: X(1)X(2)>..>X(N) *" 95/02/01 Ver.3 一、二階微分値の計算ができる(NDIF で指定)。 *" [注意] SPLINT の引数は Ver.1, Ver.2 と違う. (入力 NDIF を追加) ************************************************************************ SUBROUTINE SPLINE I ( X,Y,N,YP1,YPN, O Y2 ) IMPLICIT REAL (A-H,O-Z) PARAMETER (NMAX=200) DIMENSION X(N),Y(N), Y2(N), U(NMAX) * IF(N.GT.NMAX) THEN WRITE(6,*) '***(SPLINE)*** NOT ENOUGH WORKING AREA SIZE NMAX' STOP ENDIF * IF (YP1 .GT. .99E30) THEN !" the lower b.c. is set either to be "natural" Y2(1)=0. U(1)=0. ELSE Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) & /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE IF (YPN .GT. .99E30) THEN !" the upper b.c. is set either to be "natural" QN=0. UN=0. ELSE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END ************************************************************************ * NDIF=0:simple interpolation, NDIF=1,2:d(YA)/d(XA),d^2(YA)/d(XA)^2 * SUBROUTINE SPLINT I ( XA,YA,Y2A,N,X, NDIF, O Y ) IMPLICIT REAL (A-H,O-Z) DIMENSION XA(N),YA(N),Y2A(N) * KLO=1 KHI=N 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF((XA(N).GT.XA(1)).EQV.(XA(K).GT.X)) THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) IF (H.EQ.0.) THEN WRITE(6,*) '***(SPLINT)*** Bad XA input.' STOP ENDIF A=(XA(KHI)-X)/H B=(X-XA(KLO))/H IF (NDIF.NE.1 .AND. NDIF.NE.2) THEN Y=A*YA(KLO)+B*YA(KHI)+ & ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. ELSE IF (NDIF.EQ.1) THEN !" 1st derivative Y=(YA(KHI)-YA(KLO))/H & - (3.*A*A-1.)/6.*H*Y2A(KLO) + (3.*B*B-1.)/6.*H*Y2A(KHI) ELSE IF (NDIF.EQ.2) THEN !" 2nd derivative Y=A*Y2A(KLO)+B*Y2A(KHI) ENDIF RETURN END ************************************************************** * 相対渦度 ---> 絶対渦度 ************************************************************** SUBROUTINE CORIOL I ( HHEAD , GDATA , O HHEADS, GDATAS, D IXDIM , IYDIM , IZDIM ) * IMPLICIT NONE #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED, NFILN #else #include "gzsize.F" #endif * INTEGER IXDIM, IYDIM, IZDIM CHARACTER HHEAD ( * ) *(*) CHARACTER HHEADS ( * ) *(*) REAL GDATA ( IXDIM, IYDIM, IZDIM ) REAL GDATAS ( IXDIM, IYDIM, IZDIM ) REAL VMISS * CHARACTER HHEADA ( NDC ) * ( NCC ) INTEGER JMAX PARAMETER ( JMAX=200 ) REAL AXIS ( JMAX ) REAL F ( JMAX ) INTEGER IEOD, J, K, I REAL PI, OMEGA * IF(IYDIM.GT.JMAX) THEN CALL MSGDMP('E','CORIOL','ARRAY IS TOO SMALL') ENDIF * CALL GTSIZE(HHEADA, JMAX) CALL GUQAXV I ( HHEAD , 2 , 'LOC' , O HHEADA, AXIS , IEOD ) * PI = ATAN ( 1. )*4. OMEGA = 2. *PI / 86400. DO 2000 J=1, IYDIM F(J) = 2.*OMEGA*SIN(AXIS(J)*PI/180.) !" コリオリパラメータ 2000 CONTINUE * CALL GHPGET( HHEAD, 'MISS', VMISS ) * DO 6100 K=1,IZDIM DO 6200 J=1,IYDIM DO 6300 I=1,IXDIM IF ( GDATA( I, J, K ).NE.VMISS ) THEN GDATAS( I, J, K ) = GDATA( I, J, K ) + F ( J ) ELSE GDATAS( I, J, K ) = VMISS ENDIF 6300 CONTINUE 6200 CONTINUE 6100 CONTINUE * RETURN END ********************************************************************* *" cp from gtdsph *" ただし地球半径で割る部分を改良した。 ********************************************************************* SUBROUTINE DSPHE I ( HH , GDATA , O HHF , GDATAF, I IDIMI , JDIMI , KDIMI ) * IMPLICIT NONE INTEGER IDIMI, JDIMI, KDIMI CHARACTER HH ( * )*(*) !" ヘッダー(入力) REAL GDATA ( IDIMI*JDIMI*KDIMI ) !" データ(入力) CHARACTER HHF ( * )*(*) !" ヘッダー(出力) REAL GDATAF( IDIMI*JDIMI*KDIMI ) !" データ(出力) * LOGICAL LDX,DX !" T: x 微分 (DY=Tの方が優先される) LOGICAL LDY,DY !" T: y 微分 INTEGER IWB,WB !"(cosφ)**WB の重みを微分前に掛ける INTEGER IWA,WA !"(cosφ)**WA の重みを微分後に掛ける LOGICAL LA,A !" T: 地球半径で割る LOGICAL LZONAL,ZONAL !" T: 一番目の軸を y 軸とみなす(波数０成分用) SAVE LDX,LDY,IWB,IWA,LA,LZONAL * REAL VMISS CHARACTER HGRAD*4 * #ifdef SYS_IBMS INCLUDE (GTSINC) INCLUDE (gpainc) !" /GPAWRK/ #else #include "gtsinc.F" #include "gpainc.F" #endif *" COMMON /GPAWK?/ WDATA , ZDATA , ZWORK , *" PNM , DPNM , TRIGS , IFAX , *" NMO , GW , COSLAT, *" QPNM , QDPNM , QGW , QSINLA *" REAL ER / 6370.E3 / * *" < 配列の次元 (incl.どこが y 軸 ?) > * INTEGER IDIMS, JDIMS, KDIMS, IMAXS, JMAXS, KMAXS, JMXHF, & LMAX, MINT, MMAX, NMAX, NMDIM, & I, J, K, IJK IF(.NOT.LZONAL) THEN IDIMS=IDIMI JDIMS=JDIMI KDIMS=KDIMI IMAXS=2*(IDIMI/2) !" IDIMI が奇数なら IMAX=IDIMI-1 と仮定する JMAXS=JDIMI KMAXS=KDIMI ELSE IDIMS=1 JDIMS=IDIMI KDIMS=JDIMI IMAXS=1 JMAXS=IDIMI KMAXS=JDIMI IF(KDIMI.GT.1) CALL MSGDMP('E','DSPHE', + '3D DATA CANNOT BE PROCESSED IF -ZONAL') ENDIF * *" < 球関数の準備 > * CALL GPASPW !" 作業配列の大きさを GPASPS に知らせる I ( IDIM,JDIM,KDIM ) * CALL GPASPS !" 三角形切断しか扱えない！ I ( IMAXS , JMAXS , KMAXS , ER , O PNM , DPNM , TRIGS , IFAX , O NMO , GW , COSLAT, O LMAX , MMAX , NMAX , MINT , O NMDIM , JMXHF , W QPNM , QDPNM , QGW , QSINLA ) * *" < 微分前後の重み掛けなど > * ** IF (LA) THEN !" 地球半径で割る !" commented out by akahori ** DO 100 J=1,JMAXS ** GW(J)=GW(J)/ER ** 100 CONTINUE ** ENDIF ** CALL GHPGET( HH,'MISS',VMISS ) IF(IWB.NE.0) THEN DO 200 K=1,KDIMS DO 200 J=1,JDIMS DO 200 I=1,IDIMS IJK = I + (J-1)*IDIMS + (K-1)*IDIMS*JDIMS IF( GDATA(IJK).NE.VMISS ) THEN GDATA(IJK) = COSLAT(J)**(IWB) * GDATA(IJK) ENDIF 200 CONTINUE ENDIF * *" < 微分フラッグ > * IF(LDY) THEN HGRAD='YGRA' ELSE IF(LDX) THEN HGRAD='XGRA' ELSE HGRAD=' ' ENDIF * *" < スペクトルに (微分も) > * IF(.NOT.LZONAL) THEN CALL SPG2W M ( WDATA , I GDATA , C PNM , NMO , TRIGS , IFAX , GW , F ' ','POSO' , D IMAXS , JMAXS , KMAXS , IDIMS , JDIMS , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) ELSE CALL SPG2WM M ( WDATA , I GDATA , C PNM , NMO , TRIGS , IFAX , GW , F ' ','POSO' , D IMAXS , JMAXS , KMAXS , IDIMS , JDIMS , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) !" in this subroutine ZDATA is dummy ENDIF * *" < 格子点値に戻す > * IF(.NOT.LZONAL) THEN IF(.NOT.LDY) THEN CALL SPW2G M ( GDATAF, I WDATA , C PNM , NMO , TRIGS , IFAX , F HGRAD , 'POS' , D IMAXS , JMAXS , KMAXS , IDIMS , JDIMS , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) ELSE CALL SPW2G M ( GDATAF, I WDATA , C DPNM , NMO , TRIGS , IFAX , F HGRAD , 'POSO' , D IMAXS , JMAXS , KMAXS , IDIMS , JDIMS , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) ENDIF ELSE IF(.NOT.LDY) THEN CALL SPW2GM M ( GDATAF, I WDATA , C PNM , NMO , TRIGS , IFAX , F HGRAD , 'POS' , D IMAXS , JMAXS , KMAXS , IDIMS , JDIMS , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) ELSE CALL SPW2GM M ( GDATAF, I WDATA , C DPNM , NMO , TRIGS , IFAX , F HGRAD , 'POSO' , D IMAXS , JMAXS , KMAXS , IDIMS , JDIMS , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) ENDIF ENDIF * IF(LDY) THEN DO 400 K=1,KDIMS !"微分後の重み掛け DO 400 J=1,JDIMS DO 400 I=1,IDIMS IJK = I + (J-1)*IDIMS + (K-1)*IDIMS*JDIMS IF( GDATA(IJK).NE.VMISS ) THEN GDATAF(IJK) = COSLAT(J)**(IWA-1) * GDATAF(IJK) ENDIF 400 CONTINUE ELSE DO 500 K=1,KDIMS !"微分後の重み掛け DO 500 J=1,JDIMS DO 500 I=1,IDIMS IJK = I + (J-1)*IDIMS + (K-1)*IDIMS*JDIMS IF( GDATA(IJK).NE.VMISS ) THEN GDATAF(IJK) = COSLAT(J)**(IWA) * GDATAF(IJK) ENDIF 500 CONTINUE ENDIF IF (LA) THEN !" 地球半径で割る !" moved by akahori DO 1000 K=1,KDIMS DO 1000 J=1,JDIMS DO 1000 I=1,IDIMS IJK = I + (J-1)*IDIMS + (K-1)*IDIMS*JDIMS IF( GDATA(IJK).NE.VMISS ) THEN GDATAF(IJK) = GDATAF(IJK)/ER ENDIF 1000 CONTINUE ENDIF RETURN *=========== ENTRY EDSPHE(DX,DY,WA,WB,A,ZONAL) LDX=DX LDY=DY IWA=WA IWB=WB LA=A LZONAL=ZONAL * RETURN END ********************************************************************* *" cp from gtu2vor ********************************************************************* SUBROUTINE U2VOR I ( HU , GDU, I HV , GDV, O HVOR , GDVOR, I IMAX , JMAX , KMAX, I IMAX2 , JMAX2 , KMAX2 ) * IMPLICIT NONE INTEGER IMAX, JMAX, KMAX INTEGER IMAX2, JMAX2, KMAX2 CHARACTER HU ( * )*(*) !" ヘッダー(入力) REAL GDU ( * ) !" データ(入力) CHARACTER HV ( * )*(*) !" ヘッダー(入力) REAL GDV ( * ) !" データ(入力) CHARACTER HVOR ( * )*(*) !" ヘッダー(出力) REAL GDVOR ( * ) !" データ(出力) * #include "gzsize.F" #include "gtsinc.F" #include "gpainc.F" * REAL UVFACT ( IDIM*JDIM ) REAL WDVOR ( NMDIMC*KDIM ) !" スペクトルデータ REAL GWDEL ( JDIM ) !" 微分用のガウス荷重 INTEGER JMXHF, LMAX, MINT, MMAX, NMAX, NMDIM, & I, J, IJ, K, IJK * *" < 球関数の準備 > * CALL GPASPW !" 作業配列の大きさを GPASPS に知らせる I ( IDIM,JDIM,KDIM ) CALL GPASPS I ( IMAX , JMAX , KMAX , ER , O PNM , DPNM , TRIGS , IFAX , O NMO , GW , COSLAT, O LMAX , MMAX , NMAX , MINT , O NMDIM , JMXHF , W QPNM , QDPNM , QGW , QSINLA ) * *" < Uファクター > * DO 3210 J = 1, JMAX DO 3200 I = 1, IMAX IJ=IMAX*(J-1)+I UVFACT ( IJ ) = COSLAT(J) 3200 CONTINUE 3210 CONTINUE * *" < 微分用のガウス荷重 > * DO 3100 J = 1, JMAX GWDEL( J ) = QGW( J ) / ER / ( 1. - QSINLA(J)**2 ) 3100 CONTINUE * *" < 1. ζ のスペクトル > * DO 1100 K = 1, KMAX DO 1100 IJ = 1, IMAX*JMAX IJK=IMAX*JMAX*(K-1)+IJ GDU ( IJK ) = GDU ( IJK ) * UVFACT( IJ ) GDV ( IJK ) = GDV ( IJK ) * UVFACT( IJ ) 1100 CONTINUE * CALL SPG2W M ( WDVOR , I GDU , C DPNM , NMO , TRIGS , IFAX , GWDEL , F 'YGRA', 'POS ' , D IMAX , JMAX , KMAX , IMAX , JMAX , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) * IF ( MMAX .GE. 1 ) THEN CALL SPG2W M ( WDVOR , I GDV , C PNM , NMO , TRIGS , IFAX , GWDEL , F 'XGRA', 'ADD ' , D IMAX , JMAX , KMAX , IMAX , JMAX , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) ENDIF * *" < 格子点値に戻す > * CALL SPW2G M ( GDVOR , I WDVOR , C PNM , NMO , TRIGS , IFAX , F ' ' , 'POS ', D IMAX , JMAX , KMAX , IMAX , JMAX , D LMAX , MMAX , NMAX , MINT , NMDIM , JMXHF , W ZDATA , ZWORK ) * RETURN END ******************************************************************* *" Potential Vorticity *" = - (g/Ps) * { - ∂ｖ/∂σ * 1/(a cosφ) ∂θ/∂λ *" + ∂ｕ/∂σ * 1/a ∂θ/∂φ *" + ( f + rotｖ ) * ∂θ/∂σ } *" rot ｖ ≡ 1/(a cosφ){∂ｖ/∂λ - ∂(cosφ ｕ)/∂λ} *" 単位は PVU = 10^{-6} m^2 K / ( kg s ) ******************************************************************* SUBROUTINE PTLVOR I ( GPS , I GUS , I GVS , I GVOR , I GTHETAS, I GTHETAX, I GTHETAY, O GPV , I IMAX, JMAX, KMAX ) IMPLICIT NONE * INTEGER IMAX, JMAX, KMAX #ifdef SYS_IBMS INCLUDE (GTSINC) INCLUDE (GZSIZE) #else #include "gtsinc.F" #include "gzsize.F" #endif COMMON /GMWORK/ MWORK REAL MWORK ( IJKDIM ) * *"[入出力データの領域確保] REAL GPS (IMAX*JMAX) REAL GUS (IMAX*JMAX, KMAX) REAL GVS (IMAX*JMAX, KMAX) REAL GVOR (IMAX*JMAX, KMAX) REAL GTHETAS (IMAX*JMAX, KMAX) REAL GTHETAX (IMAX*JMAX, KMAX) REAL GTHETAY (IMAX*JMAX, KMAX) REAL GPV (IMAX*JMAX, KMAX) * REAL PVX, PVY, PVZ REAL G DATA G / 9.8 / * INTEGER IJ, K DO 10 K = 1, KMAX DO 10 IJ = 1, IMAX*JMAX PVX = - GVS(IJ,K) * GTHETAX(IJ,K) PVY = GUS(IJ,K) * GTHETAY(IJ,K) PVZ = GVOR(IJ,K)* GTHETAS(IJ,K) GPV(IJ,K) = - G / ( GPS(IJ) * 1E2 ) * ( PVX+PVY+PVZ ) * 1E6 10 CONTINUE * RETURN END