C tdraft.for C----------------------------------------------------------------------------- C 4 May 2009 - replaced REAL*8 with DOUBLE PRECISION C 13 Mar 2008 - Corrected length in common block /PCONST/ C 9 Aug 2006: Changed name F to Ftemp. Function F already exists in S8. C 28 Jul 2006: Added entry points iptempg and ptempg C 24 Apr 2006: changed 1 to 1.d0 C------------------------------------------------------------------------------ SUBROUTINE TDRAFT 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. 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 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)) 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 IPTEMPg(modtg) ! 28 Jul 2006 C------------------------------------------------------------------------------ C This subroutine is called once for each new W array. IF(W(225).NE.3.) THEN C IF(IDINT(W(225)).NE.3) THEN PRINT*,' Model check number mismatch: ' STOP ' W(225) should be 3 if you are using model TDRAFT.' ENDIF C modtg%pert_name = 'TDRAFT' 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 RETURN C C------------------------------------------------------------------------------ ENTRY IPTEMP C------------------------------------------------------------------------------ C This subroutine is called once for each new W array. IF(W(225).NE.3.) THEN C IF(IDINT(W(225)).NE.3) THEN PRINT*,' Model check number mismatch: ' STOP ' W(225) should be 3 if you are using model TDRAFT.' ENDIF C MODT(3)=8HTDRAFT MODT(4)=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 RETURN C C------------------------------------------------------------------------------ ENTRY PTEMP C------------------------------------------------------------------------------ V=0.d0 PVT=0.d0 PVR=0.d0 PVTH=0.d0 PVPH=0.d0 VR=0.d0 PVRT=0.d0 PVRR=0.d0 PVRTH=0.d0 PVRPH=0.d0 VTH=0.d0 PVTHT=0.d0 PVTHR=0.d0 PVTHTH=0.d0 PVTHPH=0.d0 VPH=0.d0 PVPHT=0.d0 PVPHR=0.d0 PVPHTH=0.d0 PVPHPH=0.d0 VH=0.d0 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) COSGAM=SINLAM0*COSTH+COSLAM0*SINTH*COSPH SINGSQ=1.d0-COSGAM**2 ! 24 Apr 2006 IF(SINGSQ.LT.0.D0) THEN WRITE (3,*) 'SINGSQ = ', SINGSQ PRINT*, 'SINGSQ', SINGSQ STOP ENDIF SINGAM=DSQRT(SINGSQ) BETA1=SINGSQ/SING1SQ BETA2=SINGSQ/SING2SQ EXBET1=DEXP(-BETA1) EXBET2=DEXP(-BETA2) SUM1=H1SQ+HSQ SUM2=H2SQ+HSQ SUM1SQ=SUM1**2 SUM2SQ=SUM2**2 TEMP=EARSQ*H*COSGAM/RSQ VR1=TEMP*(W1*EXBET1/SUM1) VR2=TEMP*(W2*EXBET2/SUM2) VR=(VR1+VR2)*H PVRR1=2.D0*VR1/R*(EARTHR*H1SQ-HCUBE)/SUM1 PVRR2=2.D0*VR2/R*(EARTHR*H2SQ-HCUBE)/SUM2 PVRR=PVRR1+PVRR2 SECGAM=1.D0/COSGAM SGPGTH=SINLAM0*SINTH-COSLAM0*COSTH*COSPH SGPGPH=COSLAM0*SINTH*SINPH PVRG1=-VR1*H*(SECGAM+2.D0*COSGAM/SING1SQ) PVRG2=-VR2*H*(SECGAM+2.D0*COSGAM/SING2SQ) PVRTH=SGPGTH*(PVRG1+PVRG2) PVRPH=SGPGPH*(PVRG1+PVRG2) TEMP1=1.D0-EXBET1 TEMP2=1.D0-EXBET2 IF(SINGAM.NE.0.D0) THEN TEMPB=-EARSQ/(R*SINGAM) VH1=TEMPB*W1*H1SQ*SING1SQ/SUM1SQ*TEMP1 VH2=TEMPB*W2*H2SQ*SING2SQ/SUM2SQ*TEMP2 VH=(VH1+VH2)*H SINALP=COSLAM0*SINPH/SINGAM SINALPSQ=SINALP**2 COSALP=(SINTH*COSTH0-SINTH0*COSTH*COSPH)/SINGAM COSALPSQ=COSALP*COSALP VTH=VH*COSALP VPH=VH*SINALP PVHR1=VH1/R*(EARTHR*H1SQ-3.D0*EARTHR*HSQ-4.D0*HCUBE)/SUM1 PVHR2=VH2/R*(EARTHR*H2SQ-3.D0*EARTHR*HSQ-4.D0*HCUBE)/SUM2 PVHR=PVHR1+PVHR2 PVTHR=PVHR*COSALP COS2ALP=1.D0-2.D0*SINALPSQ PSINALPTH=-SINALP*COSGAM/SINGSQ*SGPGTH PCOSALPTH=(COSTH*COSTH0+SINTH*SINTH0*COSPH)/SINGAM 1 -COSALP*COSGAM*SGPGTH/SINGSQ PSINALPPH=COSLAM0*COSPH/SINGAM-SINALP*COSGAM/SINGSQ*SGPGPH PCOSALPPH=SINTH0*COSTH*SINPH/SINGAM 1 -COSALP*COSGAM*SGPGPH/SINGSQ PVHPG1=VH1*H*(-1.D0+2.D0*BETA1*EXBET1/TEMP1)*COSGAM/SINGSQ PVHPG2=VH2*H*(-1.D0+2.D0*BETA2*EXBET2/TEMP2)*COSGAM/SINGSQ PVHPG=PVHPG1+PVHPG2 PVHTH=PVHPG*SGPGTH PVHPH=PVHPG*SGPGPH PVTHTH=PVHTH*COSALP+VH*PCOSALPTH PVTHPH=PVHPH*COSALP+VH*PCOSALPPH PVPHR=PVHR*SINALP PVPHTH=PVHTH*SINALP+VH*PSINALPTH PVPHPH=PVHPH*SINALP+VH*PSINALPPH ENDIF PSI1=EARSQ*W1*HSQ*SING1SQ/SUM1*TEMP1/2.D0-C1 PSI2=EARSQ*W2*HSQ*SING2SQ/SUM2*TEMP2/2.D0-C2 PSI=A1*Ftemp(A,B,PSI1,FPRIME1)+A2*Ftemp(A,B,PSI2,FPRIME2) T=T+PSI CSCGPPSI1G=R**2*H*VR1 CSCGPPSI2G=R**2*H*VR2 PPSI1R=H*W1*H1SQ*SING1SQ*TEMP1*EARSQ/SUM1SQ PPSI2R=H*W2*H2SQ*SING2SQ*TEMP2*EARSQ/SUM2SQ PPSI1TH=CSCGPPSI1G*SGPGTH PPSI2TH=CSCGPPSI2G*SGPGTH PPSI1PH=CSCGPPSI1G*SGPGPH PPSI2PH=CSCGPPSI2G*SGPGPH PTR=PTR+A1*FPRIME1*PPSI1R+A2*FPRIME2*PPSI2R PTTH=PTTH+A1*FPRIME1*PPSI1TH+A2*FPRIME2*PPSI2TH PTPH=PTPH+A1*FPRIME1*PPSI1PH+A2*FPRIME2*PPSI2PH C RETURN C------------------------------------------------------------------------------ ENTRY PTEMPg(temperature) ! 28 Jul 2006 C------------------------------------------------------------------------------ C Still need to add above formulas plus higher derivatives for gravity version RETURN END SUBROUTINE TDRAFT 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