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