C tplume.for C ----------------------------------------------------------------------!------- ! 26 Oct 2010: Added alternative definition for value of F ! 26 Oct 2010: Added common block WW to subroutine FUNC ! 26 Oct 2010: Put the calculation of F in a separate subroutine ! 25 Oct 2010: Generalized functional dependence on gamma ! 25 Oct 2010: Correctly initialized temperature perturbation and derivatives ! 9 Sep 2010: Increased to 2 plumes ! 9 Sep 2010: Moved the plume calculation to a separate subroutine ! in anticipation of adding multiple plumes C 10 Apr 2009: Compiled C 26 Feb 2009: Copied stuff from VTANH.FOR C 24 Feb 2009: Copied from tdraft2.for C ----------------------------------------------------------------------!------- SUBROUTINE TPLUME C ----------------------------------------------------------------------!------- C TEMPERATURE PERTURBATION MODEL C TEMPERATURE PERTURBATION OF A PLUME C LOCALIZED AT A SPECIFIED LONGITUDE AND LATITUDE C ----------------------------------------------------------------------!------- implicit double precision (a-h,o-z) TYPE :: values_and_derivs_1_2 DOUBLE PRECISION :: 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 DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3 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 C COMMON DECK "TT" INSERTED HERE DOUBLE PRECISION MODT COMMON/TT/MODT(4),T,PTT,PTR,PTTH,PTPH C COMMON DECK "WW" INSERTED HERE 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 (LAMBDA0,W(228)),(PH0,W(229)) & ,(LAMBDA1,W(230)),(PH1,W(231)) ! 9 Sep 2010 DOUBLE PRECISION LAMBDA0,LAMBDA1 ! 9 Sep 2010 C DOUBLE PRECISION LATITUDE,LONGITUDE C C COMMON DECK "B6" INSERTED HERE INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(10) COMMON/B6/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(611) EQUIVALENCE (DTGP,IDSDT),(ANDT,DTGP(11)) C.......The changes from the Cyber version indicated in lower case C.......below enlarge the data arrays to accommodate 100 profile points. DOUBLE PRECISION 1 ca(100),dTma(99),za(99),dla(99),DLCA(100),ALCA(100),DCA(100) 2 , cB(100),dTmB(99),zB(99),dlB(99),DLCB(100),ALCB(100),DCB(100) EQUIVALENCE (Z0A,dTGP(12)),(CA,dTgp(112)) EQUIVALENCE (Z0B,dTGP(312)),(CB,dTgp(412)) EQUIVALENCE (ZA,dTGP(13)),(dTMA,dTgp(113)),(DLA,dTgp(213)) EQUIVALENCE (ZB,dTGP(313)),(dTMB,dTgp(413)),(DLB,dTgp(513)) INTEGER dTPX ,dTQTBL(10),dTLTBL(10),dTIRMTB(10) DATA dTPX/2/ DATA dTQTBL/1,11,612,7*0/ DATA dTLTBL/1,100,8*0/ DATA dTIRMTB/1,2,8*0/ DATA RECOGdT,N/10.d0,0/ DATA AQdT/0.d0/ C ----------------------------------------------------------------------!------- ENTRY SETPTMP C ----------------------------------------------------------------------!------- ANdT=AQdT dTMX=dTPX CALL IMOVE(dTNTBL,dTQTBL,10) CALL IMOVE(dTITBL,dTLTBL,10) CALL IMOVE(dTFRMTB,dTIRMTB,10) RETURN C ----------------------------------------------------------------------!------- ENTRY IPTEMP C ----------------------------------------------------------------------!------- IF(W(225).NE.5.) THEN C IF(IDINT(W(225)).NE.5) THEN PRINT*,' Model check number mismatch: ' STOP ' W(225) should be 5 if you are using model TPLUME.' ENDIF C MODT(3)='TPLUME' MODT(4)=W(227) EARSQ=EARTHR**2 ! print *, ' lambda0 = ', lambda0 ! print *, ' ph0 = ', ph0 TH0 = PID2 - LAMBDA0 SINLAM0=DSIN(LAMBDA0) COSTH0=SINLAM0 COSLAM0=DCOS(LAMBDA0) SINTH0=COSLAM0 ! print *, ' lambda1 = ', lambda1 ! print *, ' ph1 = ', ph1 TH1 = PID2 - LAMBDA1 ! 9 Sep 2010 SINLAM1=DSIN(LAMBDA1) COSTH1=SINLAM1 COSLAM1=DCOS(LAMBDA1) SINTH1=COSLAM1 C IF HAD PREVIOUS CALL BUT NOTHING THIS TIME, EXIT NOW C RETAINING PREVIOUS TABULAR DATA COUNT IF(N.GT.0 .AND. ANdT.EQ.0.d0) RETURN C N=(ANdT+1)/6 - 2 avl=float(n) IF(N .LE. 0) 1 CALL RERROR('TPLUME','BAD N VALUE',avl) C ANdT=0.d0 C C CONVERT 'dTMA' ARRAY INPUT(OVERLAYS 'CA' ARRAY) TO 'CA' ARRAY dT0A=CA(1) CALL ITANH(N,Z0A,ZA,dT0A,dTMA,DLA,CA,CAVA,DLCA,ALCA,DCA) C C CONVERT 'dTMB' ARRAY INPUT(OVERLAYS 'CB' ARRAY) TO 'CB' ARRAY dT0B=CB(1) CALL ITANH(N,Z0B,ZB,dT0B,dTMB,DLB,CB,CAVB,DLCB,ALCB,DCB) C RETURN C C ----------------------------------------------------------------------!------- ENTRY IPTEMPg(modtg) C ----------------------------------------------------------------------!------- C This subroutine is called once for each new W array. C IF(W(225).NE.5.) THEN C IF(IDINT(W(225)).NE.5) THEN PRINT*,' Model check number mismatch: ' STOP ' W(225) should be 5 if you are using model TPLUME.' ENDIF C modtg%pert_name = 'TPLUME' modtg%pert_id = W(227) EARSQ=EARTHR**2 TH0 = PID2 - LAMBDA0 SINLAM0=DSIN(LAMBDA0) COSTH0=SINLAM0 COSLAM0=DCOS(LAMBDA0) SINTH0=COSLAM0 TH1 = PID2 - LAMBDA1 ! 9 Sep 2010 SINLAM1=DSIN(LAMBDA1) COSTH1=SINLAM1 COSLAM1=DCOS(LAMBDA1) SINTH1=COSLAM1 c NEED TO ADD SOME CALLS TO ITANH2 HERE RETURN C ----------------------------------------------------------------------!------- ENTRY PTEMP C ----------------------------------------------------------------------!------- H=R-EARTHR CALL plume(deltaT,delPTR,delPTTH,delPTPH, & H,TH,PH-PH0,TH0,SINTH0,COSLAM0,COSTH0,SINLAM0, & N,Z0A,ZA,DT0A,DLA,CA,CAVA,DLCA,ALCA,DCA, & Z0B,ZB,DT0B,DLB,CB,CAVB,DLCB,ALCB,DCB) ! 9 Sep 2010 T=T+deltaT PTR =PTR + delPTR PTTH=PTTH + delPTTH PTPH=PTPH + delPTPH if(ph1.eq.ph0.and.th1.eq.th0) return ! 9 Sep 2010 ! print *, ' lambda0 = ', lambda0 ! print *, ' ph0 = ', ph0 ! print *, ' lambda1 = ', lambda1 ! print *, ' ph1 = ', ph1 CALL plume(deltaT,delPTR,delPTTH,delPTPH, & H,TH,PH-PH1,TH1,SINTH1,COSLAM1,COSTH1,SINLAM1, & N,Z0A,ZA,DT0A,DLA,CA,CAVA,DLCA,ALCA,DCA, & Z0B,ZB,DT0B,DLB,CB,CAVB,DLCB,ALCB,DCB) ! 9 Sep 2010 T=T+deltaT PTR =PTR + delPTR PTTH=PTTH + delPTTH PTPH=PTPH + delPTPH RETURN C ----------------------------------------------------------------------!------- ENTRY PTEMPg(temperature) C ----------------------------------------------------------------------!------- C Still need to add above formulas plus higher derivatives for gravity version RETURN END SUBROUTINE TPLUME C ----------------------------------------------------------------------!------- C ----------------------------------------------------------------------!------- SUBROUTINE plume(deltaT,delPTR,delPTTH,delPTPH, & H,TH,PH,TH0,SINTH0,COSLAM0,COSTH0,SINLAM0, & N,Z0A,ZA,DT0A,DLA,CA,CAVA,DLCA,ALCA,DCA, & Z0B,ZB,DT0B,DLB,CB,CAVB,DLCB,ALCB,DCB) ! 9 Sep 2010 C ----------------------------------------------------------------------!------- implicit double precision (a-h,o-z) DOUBLE PRECISION 1 ca(100),dTma(99),za(99),dla(99),DLCA(100),ALCA(100),DCA(100) 2 , cB(100),dTmB(99),zB(99),dlB(99),DLCB(100),ALCB(100),DCB(100) COSTH=DCOS(TH) SINTH=DSIN(TH) COSPH=DCOS(PH) SINPH=DSIN(PH) C COSGAM=SINLAM0*COSTH+COSLAM0*SINTH*COSPH call bettercea(th0,th,gam,costh0,costh,cosgam, 1 sinth0,sinth,singam,ph) SINGSQ=SINGAM**2 SGPGTH=SINLAM0*SINTH-COSLAM0*COSTH*COSPH SGPGPH=COSLAM0*SINTH*SINPH deltaT=0.d0 ! 25 Oct 2010 delPTR=0.d0 ! 25 Oct 2010 delPTTH=0.d0 ! 25 Oct 2010 delPTPH=0.d0 ! 25 Oct 2010 C CALULATE STRENGTH OF TEMPERATURE PERTURBATION CALL FTANH(H,delT,PdelTR,N,Z0A,ZA,DT0A,DLA,CA,CAVA,DLCA,ALCA,DCA) if(delT.eq.0.d0) return C C CALCULATE WIDTH OF TEMPERATURE PERTURBATION CALL FTANH(H,WID,PWIDR,N,Z0B,ZB,DT0B,DLB,CB,CAVB,DLCB,ALCB,DCB) if(wid.eq.0.d0) return x=gam/wid ! 25 Oct 2010 CALL FUNC(x,F,dFdx) ! 26 Oct 2010 deltaT=delT*F pxpr = -gam*pwidr/wid**2 ! 25 Oct 2010 pFpr = dFdx*pxpr ! 25 Oct 2010 delPTR = F*pdelTr +delT*pFpr ! 25 Oct 2010 IF(SINGAM.EQ.0.) RETURN pxpCOSgam = -1.d0/(wid*singam) ! 25 Oct 2010 pFpCOSgam = dFdx*pxpCOSgam ! 25 Oct 2010 delPTTH = -SGPGTH*delT*pFpCOSgam ! 25 Oct 2010 delPTPH = -SGPGPH*delT*pFpCOSgam ! 25 Oct 2010 C RETURN C ----------------------------------------------------------------------!------- END SUBROUTINE PLUME C ----------------------------------------------------------------------!------- C ----------------------------------------------------------------------!------- SUBROUTINE FUNC(x,F,dFdx) ! 26 Oct 2010 C ----------------------------------------------------------------------!------- implicit double precision (a-h,o-z) C COMMON DECK "WW" INSERTED HERE COMMON/WW/ dum(10),MAXW,W(1000) EQUIVALENCE (EARTHR,W(1)) C DELTA TEMPERATURE 225-249 EQUIVALENCE (W(232),a2),(W(233),a4),(W(234),a6),(W(235),a8) ! 26 Oct 2010 C IF(a2.eq.0.d0.and.a4.eq.0.d0.and.a6.eq.0.d0.and.a8.eq.0.d0) then F = dexp(-x**2) ! 25 Oct 2010 dFdx = -2.d0*x*F ! 25 Oct 2010 ELSE F = dexp(-x**2*(a2+x**2*(a4+x**2*(a6+a8*x**2)))) ! 26 Oct 2010 dFdx=F*(-2.d0*x*(a2+x**2*(2.d0*a4+x**2*(3.d0*a6+4.d0*a8*x**2)))) ENDIF ! 26 Oct 2010 RETURN C ----------------------------------------------------------------------!------- END SUBROUTINE FUNC C ----------------------------------------------------------------------!-------