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