- Timestamp:
- 06/10/16 20:36:37 (8 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 1 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/etat0_dcmip2016_cyclone.f90
r388 r414 3 3 IMPLICIT NONE 4 4 PRIVATE 5 REAL(rstd),PARAMETER :: zt=15000 6 REAL(rstd),PARAMETER :: q0=0.021 7 REAL(rstd),PARAMETER :: qt=1e-11 8 REAL(rstd),PARAMETER :: T0=302.15 9 REAL(rstd),PARAMETER :: Ts=302.15 10 REAL(rstd),PARAMETER :: zq1=3000 11 REAL(rstd),PARAMETER :: zq2=8000 12 REAL(rstd),PARAMETER :: Gamma=0.007 13 REAL(rstd),PARAMETER :: pb=101500 14 REAL(rstd),PARAMETER :: Deltap=1115 15 REAL(rstd),PARAMETER :: rp=282000 16 REAL(rstd),PARAMETER :: zp=7000 17 REAL(rstd),PARAMETER :: epsilon=1e-25 18 REAL(rstd),PARAMETER :: Rd=287 19 20 REAL(rstd), PARAMETER :: lonc=Pi, latc=Pi/18, & 21 Tv0=T0*(1+0.608*q0), Tvt=Tv0-Gamma*zt 5 22 6 23 7 PUBLIC getin_etat0, compute_etat0 … … 28 12 USE mpipara, ONLY : is_mpi_root 29 13 30 IF(nqtot< 1) THEN14 IF(nqtot<5) THEN 31 15 IF (is_mpi_root) THEN 32 PRINT *, "nqtot must be at least 1 for test case DCMIP5"16 PRINT *, "nqtot must be at least 5 for test case DCMIP2016 tropical cyclone" 33 17 END IF 34 18 STOP … … 38 22 39 23 SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, Temp, ulon, ulat, q) 40 USE disvert_mod 41 USE omp_para 24 USE disvert_mod 25 USE omp_para 26 USE dcmip2016_cyclone_mod, ONLY : tropical_cyclone_test 27 USE terminator, ONLY: initial_value_Terminator 28 IMPLICIT NONE 42 29 INTEGER, INTENT(IN) :: ngrid 43 30 REAL(rstd),INTENT(IN) :: lon(ngrid) … … 50 37 REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) 51 38 39 INTEGER :: zcoords 40 REAL :: p,z 41 REAL :: rho, thetav 52 42 INTEGER :: l,ij 53 REAL(rstd) :: r, d, d1,d2, zz, Tave, Tp, vt, aa, bb54 43 55 phis(:)=0. 56 44 57 45 DO ij=1,ngrid 58 r=radius*arc(lonc,latc, lon(ij),lat(ij)) 59 ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5) 60 END DO 61 62 DO l=ll_begin,ll_end 63 aa=.5*(ap(l)+ap(l+1)) 64 bb=.5*(bp(l)+bp(l+1)) 65 DO ij=1,ngrid 66 r=radius*arc(lonc,latc, lon(ij),lat(ij)) 67 CALL compute_z(aa+bb*ps(ij),ps(ij),r,zz) 68 IF (zz<=zt) THEN 69 q(ij,l,1)=q0*exp(-zz/zq1)*exp(-(zz/zq2)**2) 70 Tave=Tv0-Gamma*zz 71 Tp=(Tv0-Gamma*zz) * ( ( 1 + 2*Rd*(Tv0-Gamma*zz)*zz / & 72 (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(zz/zp)**2))))**(-1) - 1) 73 vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & 74 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & 75 - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) 76 ELSE 77 q(ij,l,1)=qt 78 Tave=Tvt 79 Tp=0 80 vt=0 81 ENDIF 82 Temp(ij,l)=Tave+Tp 83 84 d1=sin(latc)*cos(lat(ij))-cos(latc)*sin(lat(ij))*cos(lon(ij)-lonc) 85 d2=cos(latc)*sin(lon(ij)-lonc) 86 d=MAX(1e-25,sqrt(d1**2+d2**2)) 87 ulon(ij,l)=vt*d1/d 88 ulat(ij,l)=vt*d2/d 46 z=0. 47 zcoords=1 48 CALL tropical_cyclone_test(lon(ij),lat(ij),p,z,zcoords,ulon(ij,1),ulat(ij,1), & 49 temp(ij,1),thetav,phis(ij),ps(ij),rho,q(ij,1,1)) 50 51 zcoords=0 52 DO l=ll_begin,ll_end 53 p=0.5*(ap(l)+ap(l+1) + (bp(l)+bp(l+1)) * ps(ij)) 54 CALL tropical_cyclone_test(lon(ij),lat(ij),p,z,zcoords,ulon(ij,l),ulat(ij,l), & 55 temp(ij,l),thetav,phis(ij),ps(ij),rho,q(ij,l,1)) 56 57 IF (physics_thermo==thermo_fake_moist) temp(ij,l)=Temp(ij,l)*(1+0.608*q(ij,l,1)) 58 q(ij,l,2)=0. 59 q(ij,l,3)=0. 60 CALL initial_value_Terminator(lat(ij),lon(ij),q(ij,l,4),q(ij,l,5)) 89 61 END DO 90 62 ENDDO 63 64 91 65 END SUBROUTINE compute_etat0 92 66 93 SUBROUTINE compute_z(pmodel,ps,r,z)94 USE icosa95 IMPLICIT NONE96 REAL(rstd),PARAMETER :: epsilon0=2e-1397 REAL(rstd),INTENT(IN) :: pmodel98 REAL(rstd),INTENT(IN) :: ps99 REAL(rstd),INTENT(IN) :: r100 REAL(rstd),INTENT(OUT) :: z101 102 REAL(rstd) :: pt, pave, pp, dpdz103 REAL(rstd) :: znp1104 REAL(rstd) :: epsilon105 REAL(rstd) :: p106 INTEGER :: n107 67 108 pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma))109 110 IF (pmodel>pt) THEN111 z=Tv0/Gamma*(1-(pmodel/ps)**(Rd*Gamma/g))112 ELSE113 z=zt+Rd*Tvt/g*log(Pt/pmodel)114 ENDIF115 116 IF (z>zt .OR. r>1000000) return117 118 epsilon=1119 n=0120 121 DO WHILE(epsilon>epsilon0 .AND. n<20)122 123 IF (z<zt) THEN124 pave=pb*((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))125 pp= -Deltap * exp(-(r/rp)**1.5-(z/zp)**2) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma))126 ELSE127 pave=pt*exp(g*(zt-z)/(Rd*Tvt))128 pp=0129 ENDIF130 p=pave+pp131 132 dpdz = 2*Deltap*z/zp**2 * exp(-(r/rp)**1.5) * exp(-(z/zp)**2) * ((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma)) &133 - g/(Rd*Tv0) * (preff - Deltap*exp(-(r/rp)**1.5) * exp(-(z/zp)**2)) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)-1)134 135 znp1 = z - ( pmodel - p ) / (-dpdz)136 137 epsilon=ABS( (znp1-z)/znp1 )138 z=znp1139 n=n+1140 ENDDO141 142 END SUBROUTINE compute_z143 68 144 69 END MODULE etat0_dcmip2016_cyclone_mod
Note: See TracChangeset
for help on using the changeset viewer.