| 1 | MODULE etat0_dcmip5_mod |
|---|
| 2 | USE icosa |
|---|
| 3 | PRIVATE |
|---|
| 4 | REAL(rstd),PARAMETER :: zt=15000 |
|---|
| 5 | REAL(rstd),PARAMETER :: q0=0.021 |
|---|
| 6 | REAL(rstd),PARAMETER :: qt=1e-11 |
|---|
| 7 | REAL(rstd),PARAMETER :: T0=302.15 |
|---|
| 8 | REAL(rstd),PARAMETER :: Ts=302.15 |
|---|
| 9 | REAL(rstd),PARAMETER :: zq1=3000 |
|---|
| 10 | REAL(rstd),PARAMETER :: zq2=8000 |
|---|
| 11 | REAL(rstd),PARAMETER :: Gamma=0.007 |
|---|
| 12 | REAL(rstd),PARAMETER :: pb=101500 |
|---|
| 13 | REAL(rstd),PARAMETER :: Deltap=1115 |
|---|
| 14 | REAL(rstd),PARAMETER :: rp=282000 |
|---|
| 15 | REAL(rstd),PARAMETER :: zp=7000 |
|---|
| 16 | REAL(rstd),PARAMETER :: epsilon=1e-25 |
|---|
| 17 | REAL(rstd),PARAMETER :: Rd=287 |
|---|
| 18 | |
|---|
| 19 | REAL(rstd), PARAMETER :: lonc=Pi, latc=Pi/18, & |
|---|
| 20 | Tv0=T0*(1+0.608*q0), Tvt=Tv0-Gamma*zt |
|---|
| 21 | |
|---|
| 22 | INTEGER, SAVE :: dcmip5_testcase |
|---|
| 23 | |
|---|
| 24 | PUBLIC getin_etat0, compute_etat0 |
|---|
| 25 | |
|---|
| 26 | CONTAINS |
|---|
| 27 | |
|---|
| 28 | SUBROUTINE getin_etat0 |
|---|
| 29 | dcmip5_testcase=1 |
|---|
| 30 | CALL getin("dcmip5_testcase",dcmip5_testcase) |
|---|
| 31 | END SUBROUTINE getin_etat0 |
|---|
| 32 | |
|---|
| 33 | SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, Temp, ulon, ulat, q) |
|---|
| 34 | USE disvert_mod |
|---|
| 35 | IMPLICIT NONE |
|---|
| 36 | INTEGER, INTENT(IN) :: ngrid |
|---|
| 37 | REAL(rstd),INTENT(IN) :: lon(ngrid) |
|---|
| 38 | REAL(rstd),INTENT(IN) :: lat(ngrid) |
|---|
| 39 | REAL(rstd),INTENT(OUT) :: ps(ngrid) |
|---|
| 40 | REAL(rstd),INTENT(OUT) :: phis(ngrid) |
|---|
| 41 | REAL(rstd),INTENT(OUT) :: Temp(ngrid,llm) |
|---|
| 42 | REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) |
|---|
| 43 | REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) |
|---|
| 44 | REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) |
|---|
| 45 | |
|---|
| 46 | INTEGER :: l,ij |
|---|
| 47 | REAL(rstd) :: r, d, d1,d2, zz, Tave, Tp, vt, aa, bb |
|---|
| 48 | |
|---|
| 49 | phis(:)=0. |
|---|
| 50 | |
|---|
| 51 | DO ij=1,ngrid |
|---|
| 52 | r=radius*arc(lonc,latc, lon(ij),lat(ij)) |
|---|
| 53 | ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5) |
|---|
| 54 | END DO |
|---|
| 55 | |
|---|
| 56 | DO l=1,llm |
|---|
| 57 | aa=.5*(ap(l)+ap(l+1)) |
|---|
| 58 | bb=.5*(bp(l)+bp(l+1)) |
|---|
| 59 | DO ij=1,ngrid |
|---|
| 60 | r=radius*arc(lonc,latc, lon(ij),lat(ij)) |
|---|
| 61 | CALL compute_z(aa+bb*ps(ij),ps(ij),r,zz) |
|---|
| 62 | IF (zz<=zt) THEN |
|---|
| 63 | q(ij,l,1)=q0*exp(-zz/zq1)*exp(-(zz/zq2)**2) |
|---|
| 64 | Tave=Tv0-Gamma*zz |
|---|
| 65 | Tp=(Tv0-Gamma*zz) * ( ( 1 + 2*Rd*(Tv0-Gamma*zz)*zz / & |
|---|
| 66 | (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(zz/zp)**2))))**(-1) - 1) |
|---|
| 67 | vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & |
|---|
| 68 | / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & |
|---|
| 69 | - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) |
|---|
| 70 | ELSE |
|---|
| 71 | q(ij,l,1)=qt |
|---|
| 72 | Tave=Tvt |
|---|
| 73 | Tp=0 |
|---|
| 74 | vt=0 |
|---|
| 75 | ENDIF |
|---|
| 76 | Temp(ij,l)=Tave+Tp |
|---|
| 77 | |
|---|
| 78 | d1=sin(latc)*cos(lat(ij))-cos(latc)*sin(lat(ij))*cos(lon(ij)-lonc) |
|---|
| 79 | d2=cos(latc)*sin(lon(ij)-lonc) |
|---|
| 80 | d=MAX(1e-25,sqrt(d1**2+d2**2)) |
|---|
| 81 | ulon(ij+u_right,l)=vt*d1/d |
|---|
| 82 | ulat(ij+u_right,l)=vt*d2/d |
|---|
| 83 | END DO |
|---|
| 84 | ENDDO |
|---|
| 85 | END SUBROUTINE compute_etat0 |
|---|
| 86 | |
|---|
| 87 | SUBROUTINE compute_z(pmodel,ps,r,z) |
|---|
| 88 | USE icosa |
|---|
| 89 | IMPLICIT NONE |
|---|
| 90 | REAL(rstd),PARAMETER :: epsilon0=2e-13 |
|---|
| 91 | REAL(rstd),INTENT(IN) :: pmodel |
|---|
| 92 | REAL(rstd),INTENT(IN) :: ps |
|---|
| 93 | REAL(rstd),INTENT(IN) :: r |
|---|
| 94 | REAL(rstd),INTENT(OUT) :: z |
|---|
| 95 | |
|---|
| 96 | REAL(rstd) :: pt, pave, pp, dpdz |
|---|
| 97 | REAL(rstd) :: znp1 |
|---|
| 98 | REAL(rstd) :: epsilon |
|---|
| 99 | REAL(rstd) :: p |
|---|
| 100 | INTEGER :: n |
|---|
| 101 | |
|---|
| 102 | pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma)) |
|---|
| 103 | |
|---|
| 104 | IF (pmodel>pt) THEN |
|---|
| 105 | z=Tv0/Gamma*(1-(pmodel/ps)**(Rd*Gamma/g)) |
|---|
| 106 | ELSE |
|---|
| 107 | z=zt+Rd*Tvt/g*log(Pt/pmodel) |
|---|
| 108 | ENDIF |
|---|
| 109 | |
|---|
| 110 | IF (z>zt .OR. r>1000000) return |
|---|
| 111 | |
|---|
| 112 | epsilon=1 |
|---|
| 113 | n=0 |
|---|
| 114 | |
|---|
| 115 | DO WHILE(epsilon>epsilon0 .AND. n<20) |
|---|
| 116 | |
|---|
| 117 | IF (z<zt) THEN |
|---|
| 118 | pave=pb*((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma)) |
|---|
| 119 | pp= -Deltap * exp(-(r/rp)**1.5-(z/zp)**2) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)) |
|---|
| 120 | ELSE |
|---|
| 121 | pave=pt*exp(g*(zt-z)/(Rd*Tvt)) |
|---|
| 122 | pp=0 |
|---|
| 123 | ENDIF |
|---|
| 124 | p=pave+pp |
|---|
| 125 | |
|---|
| 126 | dpdz = 2*Deltap*z/zp**2 * exp(-(r/rp)**1.5) * exp(-(z/zp)**2) * ((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma)) & |
|---|
| 127 | - g/(Rd*Tv0) * (preff - Deltap*exp(-(r/rp)**1.5) * exp(-(z/zp)**2)) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)-1) |
|---|
| 128 | |
|---|
| 129 | znp1 = z - ( pmodel - p ) / (-dpdz) |
|---|
| 130 | |
|---|
| 131 | epsilon=ABS( (znp1-z)/znp1 ) |
|---|
| 132 | z=znp1 |
|---|
| 133 | n=n+1 |
|---|
| 134 | ENDDO |
|---|
| 135 | |
|---|
| 136 | END SUBROUTINE compute_z |
|---|
| 137 | |
|---|
| 138 | END MODULE etat0_dcmip5_mod |
|---|
