C tdraft2.for C----------------------------------------------------------------------------- C 4 May 2009 - replaced REAL*8 with DOUBLE PRECISION C 13 Mar 2008 - Corrected length in common block /PCONST/ C 10 Nov 2006: Changing the formula to insure no perturbation above capping ht. C 9 Aug 2006: Correction C 8 Aug 2006: Added calculation of TH0 C 8 Aug 2006: Changed name F to Ftemp. Function F already exists in S8. C 8 Aug 2006: Correctly did the stream function in spherical polar coordinates C 4 Aug 2006: Using better formula for central earth angle C 4 Aug 2006: Formulas for small gamma (GAM) and small beta C 20 Jul 2006: Added entry points iptempg and ptempg C 10 Jul 2006: Corrected $ to % (4 places) C 27 Apr 2006: vsi -> psi C 24 Apr 2006: Modular version broken into subroutines C 24 Apr 2006: Version for the updraft being capped at a user-specified height C 24 Apr 2006: Copied from tdraft.for C------------------------------------------------------------------------------ SUBROUTINE TDRAFT2 C------------------------------------------------------------------------------ C TEMPERATURE PERTURBATION MODEL C REPRESENTS THE TEMPERATURE PERTURBATION AS A FUNCTION OF THE C STREAM FUNCTION THAT GIVES AN UPDRAFT PLUS A DOWNDRAFT. C----------------------------------------------------------------------------- implicit double precision (a-h,o-z) EXTERNAL Ftemp TYPE :: values_and_derivs_1_2 REAL(KIND=8) :: value REAL(KIND=8) :: ptm,pr,pth,pph REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm REAL(KIND=8) :: ptmr,prr,pthr,pphr REAL(KIND=8) :: ptmth,prth,pthth,pphth REAL(KIND=8) :: ptmph,prph,pthph,pphph END TYPE TYPE (values_and_derivs_1_2) :: temperature TYPE :: model CHARACTER :: model_name REAL(KIND=8) :: data_id CHARACTER :: pert_name REAL(KIND=8) :: pert_id END TYPE model TYPE (model) :: modtg REAL(KIND=8) :: h, z1, delta1, z2, delta2, c1, c2 TYPE :: value_and_derivs_1_dim REAL(KIND=8) :: value REAL(KIND=8) :: p,pp END TYPE TYPE (value_and_derivs_1_dim) :: psiA1, psiB1sg, psiC1 TYPE (value_and_derivs_1_dim) :: psiA2, psiB2sg, psiC2 TYPE :: value_and_derivs_2_dim REAL(KIND=8) :: value REAL(KIND=8) :: h,g REAL(KIND=8) :: hh,hg,gg C g means d/d cos gamma END TYPE TYPE (value_and_derivs_2_dim) :: psi1sg, psi2sg, psisg COMMON/PCONST/CREF,RGAS,GAMMA 1 ,MO2,MN2,MO,MN,MH2O,MO3 ! 9 Mar 2007 ! 13 Mar 2008 DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3 ! 13 Mar 2008 COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10 C COMMON DECK "RKAM" INSERTED HERE DOUBLE PRECISION KR,KTH,KPH COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20) C COMMON DECK "B6" INSERTED HERE INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(10) COMMON/B6/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(10) EQUIVALENCE (DTGP,IDSDT) C COMMON DECK "TT" INSERTED HERE DOUBLE PRECISION MODT COMMON/TT/MODT(4),T,PTT,PTR,PTTH,PTPH COMMON/WW/ dum(10),MAXW,W(1000) EQUIVALENCE (EARTHR,W(1)) C DELTA TEMPERATURE 225-249 EQUIVALENCE (W(225),DTMODEL),(W(226),DTFORM),(W(227),DTID) EQUIVALENCE (A1,W(228)), (A2,W(229)),(A,W(230)),(B,W(231)) EQUIVALENCE (LAMBDA0,W(128)),(PHI0,W(129)) EQUIVALENCE (W1,W(130)),(GAMMA1,W(131)),(H1,W(132)) EQUIVALENCE (W2,W(133)),(GAMMA2,W(134)),(H2,W(135)) EQUIVALENCE (z1,W(136)),(delta1,W(137)) EQUIVALENCE (z2,W(138)),(delta2,W(139)) DOUBLE PRECISION LAMBDA0 C DOUBLE PRECISION LATITUDE,LONGITUDE C COMMON DECK "B2" INSERTED HERE INTEGER DUMX,DUNTBL,DUITBL,DUFRMTBL,IDSDU(20) INTEGER DFMX,DFNTBL(10),DFITBL(10),DFFRMTB(10) DATA DFMX/1/ DATA DFNTBL/1,11,8*0/ DATA DFITBL/1,9*0/ DATA DFFRMTB/1,9*0/ C C------------------------------------------------------------------------------ ENTRY SETPTMP C------------------------------------------------------------------------------ DTMX=DFMX CALL IMOVE(DTNTBL,DFNTBL,10) CALL IMOVE(DTITBL,DFITBL,10) CALL IMOVE(DTFRMTB,DFFRMTB,10) RETURN C C------------------------------------------------------------------------------ ENTRY IPTEMP C------------------------------------------------------------------------------ IF(W(225).NE.4.) THEN C IF(IDINT(W(225)).NE.4) THEN PRINT*,' Model check number mismatch: ' STOP ' W(225) should be 4 if you are using model TDRAFT2.' ENDIF C MODT(3)=8HTDRAFT2 MODT(4)=W(227) EARSQ=EARTHR**2 SING1SQ=(DSIN(GAMMA1))**2 SING2SQ=(DSIN(GAMMA2))**2 TH0 = PID2 - LAMBDA0 ! 8 Aug 2006 SINLAM0=DSIN(LAMBDA0) COSTH0=SINLAM0 COSLAM0=DCOS(LAMBDA0) SINTH0=COSLAM0 H1SQ=H1**2 H2SQ=H2**2 C1=EARTHR*W1*SING1SQ/2.D0 ! 9 Aug 2006 C2=EARTHR*W2*SING2SQ/2.D0 ! 9 Aug 2006 CONST1 = earthr*w1/2.d0 ! 8 Aug 2006 CONST2 = earthr*w2/2.d0 ! 8 Aug 2006 RETURN C C------------------------------------------------------------------------------ ENTRY IPTEMPg(modtg) ! 20 Jul 2006 C------------------------------------------------------------------------------ C This subroutine is called once for each new W array. C IF(W(225).NE.4.) THEN C IF(IDINT(W(225)).NE.4) THEN PRINT*,' Model check number mismatch: ' STOP ' W(225) should be 4 if you are using model TDRAFT2.' ENDIF C modtg%pert_name = 'TDRAFT2' modtg%pert_id = W(227) EARSQ=EARTHR**2 SING1SQ=(DSIN(GAMMA1))**2 SING2SQ=(DSIN(GAMMA2))**2 SINLAM0=DSIN(LAMBDA0) COSTH0=SINLAM0 COSLAM0=DCOS(LAMBDA0) SINTH0=COSLAM0 H1SQ=H1**2 H2SQ=H2**2 C1=EARSQ*W1*SING1SQ/2.D0 C2=EARSQ*W2*SING2SQ/2.D0 CONST1 = earsq*w1/2.d0 CONST2 = earsq*w2/2.d0 RETURN C------------------------------------------------------------------------------ ENTRY PTEMP C------------------------------------------------------------------------------ RSQ=R**2 H=R-EARTHR HSQ=H**2 HCUBE=H*HSQ COSTH=DCOS(TH) SINTH=DSIN(TH) COSPH=DCOS(PH-PHI0) SINPH=DSIN(PH-PHI0) C COSGAM=SINLAM0*COSTH+COSLAM0*SINTH*COSPH call bettercea(th0,th,gam,costh0,costh,cosgam, 1 sinth0,sinth,singam,ph-phi0) ! 4 Aug 2006 CALL psiA(h1, h, psiA1) CALL psiB(CONST1, sing1sq, cosgam, singam, psiB1sg) ! 4 Aug 2006 CALL psiC(z1, delta1, h, psiC1) CALL psiA(h2, h, psiA2) CALL psiB(CONST2, sing2sq, cosgam, singam, psiB2sg) ! 4 Aug 2006 CALL psiC(z2, delta2, h, psiC2) psi1sg%value = (psiA1%value*psiB1sg%value - c1)*psiC1%value ! 10 Nov 2006 psi1sg%h = psiB1sg%value*(psiA1%p*psiC1%value+psiA1%value*psiC1%p) 1 -c1*psiC1%p ! 10 Nov 2006 psi1sg%g = psiA1%value*psiB1sg%p*psiC1%value c psi1sg%hh = psiB1sg%value*(psiA1%pp*psiC1%value + c 1 2.d0*psiA1%p*psiC1%p + psiA1%value*psiC1%pp) c 2 -c1*psiC1%pp ! 10 Nov 2006 c psi1sg%gg = psiA1%value*psiB1sg%pp*psiC1%value c psi1sg%hg = psiB1sg%p*(psiA1%p*psiC1%value+psiA1%value*psiC1%p) psi2sg%value = (psiA2%value*psiB2sg%value - c2)*psiC2%value ! 10 Nov 2006 psi2sg%h = psiB2sg%value*(psiA2%p*psiC2%value+psiA2%value*psiC2%p) 1 -c2*psiC2%p ! 10 Nov 2006 psi2sg%g = psiA2%value*psiB2sg%p*psiC2%value c psi2sg%hh = psiB2sg%value*(psiA2%pp*psiC2%value + c 1 2.d0*psiA2%p*psiC2%p + psiA2%value*psiC2%pp) c 2 -c2*psiC2%pp ! 10 Nov 2006 c psi2sg%gg = psiA2%value*psiB2sg%pp*psiC2%value c psi2sg%hg = psiB2sg%p*(psiA2%p*psiC2%value+psiA2%value*psiC2%p) SINGSQ=SINGAM**2 ! 4 Aug 2006 psisg%value = psi1sg%value + psi2sg%value psisg%h = psi1sg%h + psi2sg%h psisg%g = psi1sg%g + psi2sg%g c psisg%hh = psi1sg%hh + psi2sg%hh c psisg%hg = psi1sg%hg + psi2sg%hg c psisg%gg = psi1sg%gg + psi2sg%gg SGPGTH=SINLAM0*SINTH-COSLAM0*COSTH*COSPH SGPGPH=COSLAM0*SINTH*SINPH deltaT=A1*Ftemp(A,B*EARTHR,PSI1SG%value,FPRIME1) + 1 A2*Ftemp(A,B*EARTHR,PSI2SG%value,FPRIME2) T=T+deltaT PPSI1TH= -psi1sg%g *SGPGTH PPSI2TH= -psi2sg%g *SGPGTH PPSI1PH= -psi1sg%g *SGPGPH PPSI2PH= -psi2sg%g *SGPGPH PTR=PTR+A1*FPRIME1*psi1sg%h+A2*FPRIME2*psi2sg%h PTTH=PTTH+A1*FPRIME1*PPSI1TH+A2*FPRIME2*PPSI2TH PTPH=PTPH+A1*FPRIME1*PPSI1PH+A2*FPRIME2*PPSI2PH C RETURN C------------------------------------------------------------------------------ ENTRY PTEMPg(temperature) ! 20 Jul 2006 C------------------------------------------------------------------------------ C Still need to add above formulas plus higher derivatives for gravity version RETURN END SUBROUTINE TDRAFT2 C------------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION Ftemp(A,B,X,FPRIME) DOUBLE PRECISION A,B,X,FPRIME Ftemp=A+B*X FPRIME=B RETURN END FUNCTION Ftemp