C vvortx4.for C----------------------------------------------------------------------------- ! 27 Apr 2010 - Still need to add second derivatives ! 6 Apr 2010 - Moved the module to math-mod.for ! 5 Apr 2010 - Added a module for adding structures, element by element ! 2 Apr 2010 - Added structures for arguments ! 2 Apr 2010 - Added entry point IPWINDg(modug) ! 2 Apr 2010 - Added entry point PWINDg(wind) C 2 Jul 2009 - working correctly C 30 Jun 2009 - copied from vdraft2.for C----------------------------------------------------------------------------- SUBROUTINE VVORTX4 C----------------------------------------------------------------------------- C WIND VELOCITY PERTURBATION MODEL C A vortex with a viscous core and a Gaussian intensity profile in the C vertical. The axis of the vortex is vertical and may be positioned C above any geographic latitude and longitude. The vortex rotates C anticlockwise looking down. The core is essentially a solid-rotating C fluid while the outside falls off as the inverse distance from the center. C----------------------------------------------------------------------------- USE v_and_g2_maths implicit double precision (a-h,o-z) DOUBLE PRECISION :: h, z1, delta1, z2, delta2, c1, c2 DOUBLE PRECISION :: th0, gam, singam TYPE :: value_and_derivs_1_dim DOUBLE PRECISION :: value DOUBLE PRECISION :: p,pp END TYPE value_and_derivs_1_dim TYPE :: value_and_derivs_2_dim DOUBLE PRECISION :: value DOUBLE PRECISION :: h,g DOUBLE PRECISION :: hh,hg,gg C ,g means d/d cos gamma END TYPE value_and_derivs_2_dim C----------------------------------------------------------------------------- TYPE :: model CHARACTER :: model_name REAL(KIND=8) :: data_id CHARACTER :: pert_name REAL(KIND=8) :: pert_id END TYPE model TYPE (model) :: modug C----------------------------------------------------------------------------- TYPE (vector_and_gradients_1_2) :: wind, wind1, wind2 TYPE (gradient2) :: VHpp C----------------------------------------------------------------------------- C COMMON DECK "CONST" INSERTED HERE 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 COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20) double precision kr,kth,kph C COMMON DECK "UU" INSERTED HERE DOUBLE PRECISION MODU, LAMBDA0 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 (u0,W(130)),(GAMMA0,W(131)),(Hmax,W(132)) EQUIVALENCE (WH,W(133)) 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 C----------------------------------------------------------------------------- ENTRY SETPWN C----------------------------------------------------------------------------- c print*,'Beginning of setpwn' DUMX=DUPX CALL IMOVE(DUNTBL,DUQTBL,10) CALL IMOVE(DUITBL,DULTBL,10) CALL IMOVE(DUFRMTBL,DUIRMTBL,10) RETURN C C----------------------------------------------------------------------------- ENTRY IPWINDR C----------------------------------------------------------------------------- IF(W(125).NE.3.) THEN PRINT*,' Model check number mismatch: ' STOP ' W(125) should be 3 if you are using model VVORTX4.' ENDIF MODU(3)=8HVVORTX4 MODU(4)=W(127) SING0SQ=(DSIN(GAMMA0))**2 th0 = pid2 - Lambda0 SINLAM0=DSIN(LAMBDA0) COSTH0=SINLAM0 COSLAM0=DCOS(LAMBDA0) SINTH0=COSLAM0 HmSQ=Hmax**2 WHSQ=WH**2 DENOM = 0.D0 IF(WHSQ.NE.0.D0) DENOM=1.D0/WHSQ C2 = 1.256431209 C1 = 1.D0+1.D0/(2.D0*C2) C3 = C1*U0*GAMMA0 C4 = C2/GAMMA0**2 RETURN C C----------------------------------------------------------------------------- ENTRY IPWINDg(modug) C----------------------------------------------------------------------------- IF(W(125).NE.3.) THEN PRINT*,' Model check number mismatch: ' STOP ' W(125) should be 3 if you are using model VVORTX4.' ENDIF MODU(3)=8HVVORTX4 MODU(4)=W(127) modug%pert_name = 'VVORTX4' modug%pert_id = W(127) SING0SQ=(DSIN(GAMMA0))**2 th0 = pid2 - Lambda0 SINLAM0=DSIN(LAMBDA0) COSTH0=SINLAM0 COSLAM0=DCOS(LAMBDA0) SINTH0=COSLAM0 HmSQ=Hmax**2 WHSQ=WH**2 DENOM = 0.D0 IF(WHSQ.NE.0.D0) DENOM=1.D0/WHSQ C2 = 1.256431209 C1 = 1.D0+1.D0/(2.D0*C2) C3 = C1*U0*GAMMA0 C4 = C2/GAMMA0**2 RETURN C C----------------------------------------------------------------------------- ENTRY PWINDR C----------------------------------------------------------------------------- H=R-EARTHR 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) GAMMASQ=GAM**2 SINGSQ=SINGAM**2 SGPGTH=SINLAM0*SINTH-COSLAM0*COSTH*COSPH SGPGPH=COSLAM0*SINTH*SINPH VR=0.D0 PVRR=0.D0 PVRTH = 0.D0 PVRPH = 0.D0 IF(SINGAM.NE.0.D0) THEN A=C3/GAM D=DEXP(-C4*GAMMASQ) B=(1.D0-D) C=DEXP(-DENOM*(H-Hmax)**2) VH=A*B*C PAPG=-A/GAM PBPG=2.D0*C4*GAM*D PVHPG=-(A*PBPG+B*PAPG)*C/SINGAM PCPR =-2.D0*(H-Hmax)*DENOM*C PVHR=A*B*PCPR SINALP=COSLAM0*SINPH/SINGAM SINALPSQ=SINALP**2 COSALP=(SINTH*COSTH0-SINTH0*COSTH*COSPH)/SINGAM c cosalp = (COSTH0-COSTH*COSGAM)/(SINTH*SINGAM) ! equivalent COSALPSQ=COSALP*COSALP c one=sinalpsq+cosalpsq VPH=VH*COSALP VTH=-VH*SINALP PVTHR=-PVHR*SINALP 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 PVHTH=-PVHPG*SGPGTH PVHPH=-PVHPG*SGPGPH PVPHR= PVHR*COSALP PVPHTH=PVHTH*COSALP+VH*PCOSALPTH PVPHPH=PVHPH*COSALP+VH*PCOSALPPH PVTHR=-PVHR*SINALP PVTHTH=-PVHTH*SINALP-VH*PSINALPTH PVTHPH=-PVHPH*SINALP-VH*PSINALPPH ELSE VTH=0.d0 VPH=0.D0 PVTHR=0.d0 PVPHR=0.d0 PVTHTH=0.D0 PVPHTH=C3*C4*C PVTHPH=-C3*C4*C*SINTH0 PVPHPH=0.D0 ENDIF C PVRT=0.d0 PVTHT=0.d0 PVPHT=0.d0 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 c print *, ' vr0 = ', vr0 c print *, ' vth0 = ', vth0 c print *, ' vph0 = ', vph0 VSQ0 = VR0**2 + VTH0**2 + VPH0**2 V0 = DSQRT(VSQ0) IF (V0 /= 0.D0) THEN 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 ELSE PVT0 = 0.d0 PVR0 = 0.d0 PVTH0 = 0.d0 PVPH0 = 0.d0 ENDIF C print *, ' pvr0 = ', pvr0 C print *, ' pvth0 = ', pvth0 C print *, ' pvph0 = ', pvph0 C RETURN C----------------------------------------------------------------------------- ENTRY PWINDg(wind) C----------------------------------------------------------------------------- H=R-EARTHR 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) GAMMASQ=GAM**2 SINGSQ=SINGAM**2 SGPGTH=SINLAM0*SINTH-COSLAM0*COSTH*COSPH SGPGPH=COSLAM0*SINTH*SINPH VR=0.D0 PVRR=0.D0 PVRTH = 0.D0 PVRPH = 0.D0 IF(SINGAM.NE.0.D0) THEN A=C3/GAM D=DEXP(-C4*GAMMASQ) B=(1.D0-D) C=DEXP(-DENOM*(H-Hmax)**2) VH=A*B*C PAPG=-A/GAM PBPG=2.D0*C4*GAM*D PVHPG=-(A*PBPG+B*PAPG)*C/SINGAM PCPR =-2.D0*(H-Hmax)*DENOM*C PVHR=A*B*PCPR SINALP=COSLAM0*SINPH/SINGAM SINALPSQ=SINALP**2 COSALP=(SINTH*COSTH0-SINTH0*COSTH*COSPH)/SINGAM c cosalp = (COSTH0-COSTH*COSGAM)/(SINTH*SINGAM) ! equivalent COSALPSQ=COSALP*COSALP c one=sinalpsq+cosalpsq VPH=VH*COSALP VTH=-VH*SINALP PVTHR=-PVHR*SINALP 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 PVHTH=-PVHPG*SGPGTH PVHPH=-PVHPG*SGPGPH PVPHR= PVHR*COSALP PVPHTH=PVHTH*COSALP+VH*PCOSALPTH PVPHPH=PVHPH*COSALP+VH*PCOSALPPH PVTHR=-PVHR*SINALP PVTHTH=-PVHTH*SINALP-VH*PSINALPTH PVTHPH=-PVHPH*SINALP-VH*PSINALPPH ELSE VTH=0.d0 VPH=0.D0 PVTHR=0.d0 PVPHR=0.d0 PVTHTH=0.D0 PVPHTH=C3*C4*C PVTHPH=-C3*C4*C*SINTH0 PVPHPH=0.D0 ENDIF C PVRT=0.d0 PVTHT=0.d0 PVPHT=0.d0 wind1%r%v = VR wind1%r%p%tm = PVRT wind1%r%p%r = PVRR wind1%r%p%th = PVRTH wind1%r%p%ph = PVRPH wind1%r%pp%tm%tm = 0.d0 wind1%r%pp%tm%r = 0.d0 wind1%r%pp%tm%th = 0.d0 wind1%r%pp%tm%ph = 0.d0 wind1%r%pp%r%tm = 0.d0 wind1%r%pp%r%r = 0.d0 wind1%r%pp%r%th = 0.d0 wind1%r%pp%r%ph = 0.d0 wind1%r%pp%th%tm = 0.d0 wind1%r%pp%th%r = 0.d0 wind1%r%pp%th%th = 0.d0 wind1%r%pp%th%ph = 0.d0 wind1%r%pp%ph%tm = 0.d0 wind1%r%pp%ph%r = 0.d0 wind1%r%pp%ph%th = 0.d0 wind1%r%pp%ph%ph = 0.d0 wind1%th%v = VTH wind1%th%p%tm = PVTHT wind1%th%p%r = PVTHR wind1%th%p%th = PVTHTH wind1%th%p%ph = PVTHPH wind1%th%pp%tm%tm = 0.d0 wind1%th%pp%tm%r = 0.d0 wind1%th%pp%tm%th = 0.d0 wind1%th%pp%tm%ph = 0.d0 wind1%th%pp%r%tm = 0.d0 wind1%th%pp%r%r = 0.d0 wind1%th%pp%r%th = 0.d0 wind1%th%pp%r%ph = 0.d0 wind1%th%pp%th%tm = 0.d0 wind1%th%pp%th%r = 0.d0 wind1%th%pp%th%th = 0.d0 wind1%th%pp%th%ph = 0.d0 wind1%th%pp%ph%tm = 0.d0 wind1%th%pp%ph%r = 0.d0 wind1%th%pp%ph%th = 0.d0 wind1%th%pp%ph%ph = 0.d0 wind1%ph%v = VPH wind1%ph%p%tm = PVPHT wind1%ph%p%r = PVPHR wind1%ph%p%th = PVPHTH wind1%ph%p%ph = PVPHPH wind1%ph%pp%tm%tm = 0.d0 wind1%ph%pp%tm%r = 0.d0 wind1%ph%pp%tm%th = 0.d0 wind1%ph%pp%tm%ph = 0.d0 wind1%ph%pp%r%tm = 0.d0 wind1%ph%pp%r%r = 0.d0 wind1%ph%pp%r%th = 0.d0 wind1%ph%pp%r%ph = 0.d0 wind1%ph%pp%th%tm = 0.d0 wind1%ph%pp%th%r = 0.d0 wind1%ph%pp%th%th = 0.d0 wind1%ph%pp%th%ph = 0.d0 wind1%ph%pp%ph%tm = 0.d0 wind1%ph%pp%ph%r = 0.d0 wind1%ph%pp%ph%th = 0.d0 wind1%ph%pp%ph%ph = 0.d0 C----------------------------------------------------------------------------- C don't need below now. wind1%A%v = V wind1%A%p%tm = PVT wind1%A%p%r = PVR wind1%A%p%th = PVTH wind1%A%p%ph = PVPH wind1%A%v = 0.d0 wind1%A%p%tm = 0.d0 wind1%A%p%r = 0.d0 wind1%A%p%th = 0.d0 wind1%A%p%ph = 0.d0 wind1%A%pp%tm%tm = 0.d0 wind1%A%pp%tm%r = 0.d0 wind1%A%pp%tm%th = 0.d0 wind1%A%pp%tm%ph = 0.d0 wind1%A%pp%r%tm = 0.d0 wind1%A%pp%r%r = 0.d0 wind1%A%pp%r%th = 0.d0 wind1%A%pp%r%ph = 0.d0 wind1%A%pp%th%tm = 0.d0 wind1%A%pp%th%r = 0.d0 wind1%A%pp%th%th = 0.d0 wind1%A%pp%th%ph = 0.d0 wind1%A%pp%ph%tm = 0.d0 wind1%A%pp%ph%r = 0.d0 wind1%A%pp%ph%th = 0.d0 wind1%A%pp%ph%ph = 0.d0 C----------------------------------------------------------------------------- wind2 = wind + wind1 wind = wind2 V0 = wind%A%v VR0 = wind%r%v PVRT0 = wind%r%p%tm PVRR0 = wind%r%p%r PVRTH0 = wind%r%p%th PVRPH0 = wind%r%p%ph VTH0 = wind%th%v PVTHT0 = wind%th%p%tm PVTHR0 = wind%th%p%r PVTHTH0 = wind%th%p%th PVTHPH0 = wind%th%p%ph VPH0 = wind%ph%v PVPHT0 = wind%ph%p%tm PVPHR0 = wind%ph%p%r PVPHTH0 = wind%ph%p%th PVPHPH0 = wind%ph%p%ph c print *, ' vr0 = ', vr0 c print *, ' vth0 = ', vth0 c print *, ' vph0 = ', vph0 C----------------------------------------------------------------------------- C don't need below now. IF (V0 /= 0.D0) THEN wind%A%p%tm = (VR0*PVRT0 + VTH0*PVTHT0 + VPH0*PVPHT0)/V0 wind%A%p%r = (VR0*PVRR0 + VTH0*PVTHR0 + VPH0*PVPHR0)/V0 wind%A%p%th = (VR0*PVRTH0 + VTH0*PVTHTH0 + VPH0*PVPHTH0)/V0 wind%A%p%ph = (VR0*PVRPH0 + VTH0*PVTHPH0 + VPH0*PVPHPH0)/V0 ELSE wind%A%p%tm = 0.d0 wind%A%p%r = 0.d0 wind%A%p%th = 0.d0 wind%A%p%ph = 0.d0 ENDIF C----------------------------------------------------------------------------- V0 = wind%A%v PVT0 = wind%A%p%tm PVR0 = wind%A%p%r PVTH0 = wind%A%p%th PVPH0 = wind%A%p%ph C print *, ' pvr0 = ', pvr0 C print *, ' pvth0 = ', pvth0 C print *, ' pvph0 = ', pvph0 RETURN C----------------------------------------------------------------------------- END SUBROUTINE VVORTX4 C-----------------------------------------------------------------------------