C vdraft.for C----------------------------------------------------------------------------- C 13 Mar 2008 - Corrected length in common block /PCONST/ C 10 Mar 2008 - corrected length errors in blank common C 24 Apr 2006: changed 1 to 1.d0 and 0 to 0.d0 C 15 Jun 2005: Converted to perturbation model C 31 May 2005: Converted back to old check for wrong model C 28 May 2005: Commented out old programming C 27 May 2005: Deleted old programming C 27 May 2005: Commented out old programming C 27 May 2005: Corrected PVHR1 and PVHR2 C 27 May 2005: Corrected SING1SQ and SING2SQ C 26 May 2005: Corrected PCOSALPTH and PCOSALPPH C 26 May 2005: Print out phi, theta, alpha C 25 May 2005: Corrected cos alpha calculation plus derivatives. C 24 May 2005: Added printing of gradients. C 23 May 2005: Reverted to new formulas. C 23 May 2005: Reverted to old alternate without errors C 23 May 2005: Reverted to cos^2 alpha = 1.d-16 when negative. C 23 May 2005: Reverted to cos^2 alpha = zero when negative. C 23 May 2005: Reverted to cos^2 alpha = 1.d-16 when negative. C 23 May 2005: Reverted to cos^2 alpha = 1.d-38 when negative. C 23 May 2005: Reverted to old alternate with errors C 19 May 2005: included both old and new formulas C 18 May 2005: Corrected formulas for PVTHTH and PVPHTH C 18 May 2005: Corrected formulas for PVTHPH and PVPHPH C 17 May 2005: Added and changed some debugging statements. C 16 May 2005: Added some debugging statements. C 14 May 2005: Corrected formula for PVHR C 12 May 2005: Reorganized to avoid dividing by zero. C 11 May 2005: Added print statements at end C 3 May 2005: Declared LAMBDA0 double precision. C----------------------------------------------------------------------------- SUBROUTINE VDRAFT C----------------------------------------------------------------------------- C WIND VELOCITY MODEL C An updraft plus a downdraft centered at a specified longitude and C latitude, based on a stream function. implicit double precision (a-h,o-z) 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 COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20) double precision kr,kth,kph ! 10 Mar 2008 DOUBLE PRECISION MODU, LAMBDA0 ! 3 May 2005 COMMON/UU/MODU(4) 1 ,V0 ,PVT0 ,PVR0 ,PVTH0 ,PVPH0 2 ,VR0 ,PVRT0 ,PVRR0 ,PVRTH0 ,PVRPH0 3 ,VTH0,PVTHT0,PVTHR0,PVTHTH0,PVTHPH0 4 ,VPH0,PVPHT0,PVPHR0,PVPHTH0,PVPHPH0 COMMON/WW/ dum(10), MAXW,W(1000) EQUIVALENCE (EARTHR,W(1)) C WIND PERTURBATION 125-149 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)) C COMMON DECK "B2" INSERTED HERE INTEGER DUMX,DUNTBL,DUITBL,DUFRMTBL,IDSDU(20) COMMON/B2/DUMX,DUNTBL(10),DUITBL(10),DUFRMTBL(10),DUGP(10) EQUIVALENCE (DUGP,IDSDU) INTEGER DUPX ,DUQTBL(10),DULTBL(10),DUIRMTBL(10) DATA DUPX/1/ DATA DUQTBL/1,11,8*0/ DATA DULTBL/1,9*0/ DATA DUIRMTBL/1,9*0/ C ENTRY SETPWN DUMX=DUPX CALL IMOVE(DUNTBL,DUQTBL,10) CALL IMOVE(DUITBL,DULTBL,10) CALL IMOVE(DUFRMTBL,DUIRMTBL,10) RETURN C ENTRY IPWINDR IF(W(125).NE.1.) THEN C IF(IDINT(W(125)).NE.1) THEN PRINT*,' Model check number mismatch: ' STOP ' W(125) should be 1 if you are using model VDRAFT.' ENDIF MODU(3)=8HVDRAFT MODU(4)=W(127) EARSQ=EARTHR**2 SING1SQ=(DSIN(GAMMA1))**2 ! 27 May 2005 SING2SQ=(DSIN(GAMMA2))**2 ! 27 May 2005 SINLAM0=DSIN(LAMBDA0) COSTH0=SINLAM0 ! 25 May 2005: COSLAM0=DCOS(LAMBDA0) SINTH0=COSLAM0 ! 25 May 2005: H1SQ=H1**2 H2SQ=H2**2 C CALL IPWINDR RETURN C ENTRY PWINDR V=0.d0 ! 24 Apr 2006 and below 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 ! 16 May 2005 stop ! 16 May 2005 ENDIF ! 16 May 2005 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) IF(SINGAM.NE.0.D0) THEN TEMPB=-EARSQ/(R*SINGAM) TEMP1=1.D0-EXBET1 TEMP2=1.D0-EXBET2 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 ! 25 May 2005 COSALPSQ=COSALP*COSALP ! 25 May 2005 one=sinalpsq+cosalpsq ! 26 May 2005 VTH=VH*COSALP VPH=VH*SINALP PVHR1=VH1/R*(EARTHR*H1SQ-3.D0*EARTHR*HSQ-4.D0*HCUBE)/SUM1 ! 27 May 2005 PVHR2=VH2/R*(EARTHR*H2SQ-3.D0*EARTHR*HSQ-4.D0*HCUBE)/SUM2 ! 27 May 2005 PVHR=PVHR1+PVHR2 PVTHR=PVHR*COSALP COS2ALP=1.D0-2.D0*SINALPSQ ! 17 May 2005 PSINALPTH=-SINALP*COSGAM/SINGSQ*SGPGTH ! 18 May 2005 PCOSALPTH=(COSTH*COSTH0+SINTH*SINTH0*COSPH)/SINGAM ! 25 May 2005 1 -COSALP*COSGAM*SGPGTH/SINGSQ ! 26 May 2005 PSINALPPH=COSLAM0*COSPH/SINGAM-SINALP*COSGAM/SINGSQ*SGPGPH !18 May 2005 PCOSALPPH=SINTH0*COSTH*SINPH/SINGAM 1 -COSALP*COSGAM*SGPGPH/SINGSQ ! 26 May 2005 PVHPG1=VH1*H*(-1.D0+2.D0*BETA1*EXBET1/TEMP1)*COSGAM/SINGSQ !18 May 2005 PVHPG2=VH2*H*(-1.D0+2.D0*BETA2*EXBET2/TEMP2)*COSGAM/SINGSQ !18 May 2005 PVHPG=PVHPG1+PVHPG2 PVHTH=PVHPG*SGPGTH PVHPH=PVHPG*SGPGPH PVTHTH=PVHTH*COSALP+VH*PCOSALPTH ! 18 May 2005 PVTHPH=PVHPH*COSALP+VH*PCOSALPPH ! 18 May 2005 PVPHR=PVHR*SINALP PVPHTH=PVHTH*SINALP+VH*PSINALPTH ! 18 May 2005 PVPHPH=PVHPH*SINALP+VH*PSINALPPH ! 18 May 2005 ENDIF C VSQ=VR**2+VH**2 C IF(VSQ.Lt.0.D0) stop C V=DSQRT(VSQ) C IF(V.NE.0.D0) THEN C PVTH=(VR*PVRTH + VTH*PVTHTH + VPH*PVPHTH)/V C PVPH=(VR*PVRPH + VTH*PVTHPH + VPH*PVPHPH)/V C PVR=(VR*PVRR + VTH*PVTHR + VPH*PVPHR)/V C ENDIF C CALL PWINDR C VR0 = VR0 + VR PVRT0 = PVRT0 + PVRT PVRR0 = PVRR0 + PVRR PVRTH0 = PVRTH0 + PVRTH PVRPH0 = PVRPH0 + PVRPH VTH0 = VTH0 + VTH PVTHT0 = PVTHT0 + PVTHT PVTHR0 = PVTHR0 + PVTHR PVTHTH0 = PVTHTH0 + PVTHTH PVTHPH0 = PVTHPH0 + PVTHPH VPH0 = VPH0 + VPH PVPHT0 = PVPHT0 + PVPHT PVPHR0 = PVPHR0 + PVPHR PVPHTH0 = PVPHTH0 + PVPHTH PVPHPH0 = PVPHPH0 + PVPHPH C VSQ0 = VR0**2 + VTH0**2 + VPH0**2 IF(VSQ0.LT.0.D0) THEN PRINT*, 'Velocity cannot be zero!' STOP ENDIF V0 = DSQRT(VSQ0) PVT0 = 0.d0 PVR0 = 0.d0 PVTH0 = 0.d0 PVPH0 = 0.d0 IF (V0.EQ.0.d0) RETURN PVT0 = (VR0*PVRT0 + VTH0*PVTHT0 + VPH0*PVPHT0)/V0 PVR0 = (VR0*PVRR0 + VTH0*PVTHR0 + VPH0*PVPHR0)/V0 PVTH0 = (VR0*PVRTH0 + VTH0*PVTHTH0 + VPH0*PVPHTH0)/V0 PVPH0 = (VR0*PVRPH0 + VTH0*PVTHPH0 + VPH0*PVPHPH0)/V0 C RETURN END SUBROUTINE VDRAFT