[113] | 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 | REAL(rstd) :: lonc |
---|
| 19 | REAL(rstd) :: latc |
---|
| 20 | TYPE(t_field),POINTER :: f_out_i(:) |
---|
| 21 | REAL(rstd),POINTER :: out_i(:,:) |
---|
| 22 | |
---|
| 23 | PUBLIC etat0 |
---|
| 24 | CONTAINS |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | |
---|
| 28 | SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) |
---|
| 29 | USE icosa |
---|
| 30 | IMPLICIT NONE |
---|
| 31 | TYPE(t_field),POINTER :: f_ps(:) |
---|
| 32 | TYPE(t_field),POINTER :: f_phis(:) |
---|
| 33 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
---|
| 34 | TYPE(t_field),POINTER :: f_u(:) |
---|
| 35 | TYPE(t_field),POINTER :: f_q(:) |
---|
| 36 | |
---|
| 37 | REAL(rstd),POINTER :: ps(:) |
---|
| 38 | REAL(rstd),POINTER :: phis(:) |
---|
| 39 | REAL(rstd),POINTER :: theta_rhodz(:,:) |
---|
| 40 | REAL(rstd),POINTER :: u(:,:) |
---|
| 41 | REAL(rstd),POINTER :: q(:,:,:) |
---|
| 42 | INTEGER :: ind |
---|
| 43 | |
---|
| 44 | CALL allocate_field(f_out_i,field_t,type_real,llm) |
---|
| 45 | |
---|
| 46 | DO ind=1,ndomain |
---|
| 47 | CALL swap_dimensions(ind) |
---|
| 48 | CALL swap_geometry(ind) |
---|
| 49 | ps=f_ps(ind) |
---|
| 50 | phis=f_phis(ind) |
---|
| 51 | theta_rhodz=f_theta_rhodz(ind) |
---|
| 52 | u=f_u(ind) |
---|
| 53 | q=f_q(ind) |
---|
| 54 | out_i=f_out_i(ind) |
---|
| 55 | CALL compute_etat0_dcmip5(ps, phis, theta_rhodz, u, q) |
---|
| 56 | ENDDO |
---|
| 57 | CALL writefield("out_i",f_out_i) |
---|
| 58 | CALL deallocate_field(f_out_i) |
---|
| 59 | END SUBROUTINE etat0 |
---|
| 60 | |
---|
| 61 | SUBROUTINE compute_etat0_dcmip5(ps, phis, theta_rhodz, ue, q) |
---|
| 62 | USE icosa |
---|
| 63 | USE pression_mod |
---|
| 64 | USE wind_mod |
---|
| 65 | USE theta2theta_rhodz_mod |
---|
[114] | 66 | USE exner_mod |
---|
[113] | 67 | IMPLICIT NONE |
---|
| 68 | REAL(rstd),INTENT(OUT) :: ps(iim*jjm) |
---|
| 69 | REAL(rstd),INTENT(OUT) :: phis(iim*jjm) |
---|
| 70 | REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) |
---|
| 71 | REAL(rstd),INTENT(OUT) :: ue(3*iim*jjm,llm) |
---|
| 72 | REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) |
---|
| 73 | |
---|
| 74 | INTEGER :: i,j,l,ij |
---|
| 75 | INTEGER :: testcase |
---|
| 76 | REAL(rstd) :: Tv0, Tvt |
---|
| 77 | REAL(rstd) :: r |
---|
| 78 | REAL(rstd) :: z(iim*jjm,llm),zz |
---|
| 79 | REAL(rstd) :: p(iim*jjm,llm+1) |
---|
| 80 | REAL(rstd) :: Tave,Tp |
---|
| 81 | REAL(rstd) :: T(iim*jjm,llm) |
---|
| 82 | REAL(rstd) :: vt |
---|
| 83 | REAL(rstd) :: d1,d2,d |
---|
| 84 | REAL(rstd) :: ulon(3*iim*jjm,llm) |
---|
| 85 | REAL(rstd) :: ulat(3*iim*jjm,llm) |
---|
| 86 | REAL(rstd) :: lon,lat |
---|
[114] | 87 | REAL(rstd) :: pks(iim*jjm) |
---|
| 88 | REAL(rstd) :: pk(iim*jjm,llm) |
---|
[113] | 89 | |
---|
| 90 | |
---|
| 91 | latc=Pi/18 |
---|
| 92 | lonc=Pi |
---|
| 93 | |
---|
| 94 | Tv0=T0/(1+0.608*q0) |
---|
| 95 | Tvt=Tv0-Gamma*zt |
---|
| 96 | |
---|
| 97 | testcase=1 |
---|
| 98 | CALL getin("dcmip5_testcase",testcase) |
---|
| 99 | |
---|
| 100 | phis(:)=0. |
---|
| 101 | |
---|
| 102 | DO j=jj_begin,jj_end |
---|
| 103 | DO i=ii_begin,ii_end |
---|
| 104 | ij=(j-1)*iim+i |
---|
| 105 | CALL xyz2lonlat(xyz_i(ij,:),lon,lat) |
---|
| 106 | r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) |
---|
| 107 | ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5) |
---|
| 108 | ENDDO |
---|
| 109 | ENDDO |
---|
| 110 | |
---|
| 111 | CALL compute_pression(ps,p,0) |
---|
| 112 | |
---|
| 113 | DO l=1,llm |
---|
| 114 | DO j=jj_begin,jj_end |
---|
| 115 | DO i=ii_begin,ii_end |
---|
| 116 | ij=(j-1)*iim+i |
---|
| 117 | CALL xyz2lonlat(xyz_i(ij,:),lon,lat) |
---|
| 118 | r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) |
---|
| 119 | CALL compute_z(0.5*(p(ij,l)+p(ij,l+1)),ps(ij),r,z(ij,l)) |
---|
| 120 | ENDDO |
---|
| 121 | ENDDO |
---|
| 122 | ENDDO |
---|
| 123 | out_i=z |
---|
| 124 | |
---|
| 125 | DO l=1,llm |
---|
| 126 | DO j=jj_begin-1,jj_end+1 |
---|
| 127 | DO i=ii_begin-1,ii_end+1 |
---|
| 128 | ij=(j-1)*iim+i |
---|
| 129 | IF (z(ij,l)<=zt) THEN |
---|
| 130 | q(ij,l,1)=q0*exp(-z(ij,l)/zq1)*exp(-(z(ij,l)/zq2)**2) |
---|
| 131 | ELSE |
---|
| 132 | q(ij,l,1)=qt |
---|
| 133 | ENDIF |
---|
| 134 | ENDDO |
---|
| 135 | ENDDO |
---|
| 136 | ENDDO |
---|
| 137 | |
---|
| 138 | |
---|
| 139 | DO l=1,llm |
---|
| 140 | DO j=jj_begin-1,jj_end+1 |
---|
| 141 | DO i=ii_begin-1,ii_end+1 |
---|
| 142 | ij=(j-1)*iim+i |
---|
| 143 | CALL xyz2lonlat(xyz_i(ij,:),lon,lat) |
---|
| 144 | r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) |
---|
| 145 | IF (z(ij,l)<=zt) THEN |
---|
| 146 | Tave=Tv0-Gamma*z(ij,l) |
---|
| 147 | Tp=(Tv0-Gamma*z(ij,l)) * ( ( 1 + 2*Rd*(Tv0-Gamma*z(ij,l))*z(ij,l) / & |
---|
| 148 | (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(z(ij,l)/zp)**2))))**(-1) - 1) |
---|
| 149 | ELSE |
---|
| 150 | Tave=Tvt |
---|
| 151 | Tp=0 |
---|
| 152 | ENDIF |
---|
| 153 | T(ij,l)=Tave+Tp |
---|
[114] | 154 | out_i(ij,l)=Tave+Tp |
---|
[113] | 155 | ENDDO |
---|
| 156 | ENDDO |
---|
| 157 | ENDDO |
---|
| 158 | |
---|
| 159 | DO l=1,llm |
---|
| 160 | DO j=jj_begin-1,jj_end+1 |
---|
| 161 | DO i=ii_begin-1,ii_end+1 |
---|
| 162 | ij=(j-1)*iim+i |
---|
| 163 | |
---|
| 164 | CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) |
---|
| 165 | zz=0.5*(z(ij,l)+z(ij+t_right,l)) |
---|
| 166 | r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) |
---|
| 167 | IF (zz<=zt) THEN |
---|
| 168 | vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & |
---|
| 169 | / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & |
---|
| 170 | - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) |
---|
| 171 | ELSE |
---|
| 172 | vt = 0 |
---|
| 173 | ENDIF |
---|
| 174 | d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc) |
---|
| 175 | d2=cos(latc)*sin(lon-lonc) |
---|
| 176 | d=MAX(1e-25,sqrt(d1**2+d2**2)) |
---|
| 177 | ulon(ij+u_right,l)=vt*d1/d |
---|
| 178 | ulat(ij+u_right,l)=vt*d2/d |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) |
---|
| 183 | zz=0.5*(z(ij,l)+z(ij+t_lup,l)) |
---|
| 184 | r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) |
---|
| 185 | IF (zz<=zt) THEN |
---|
| 186 | vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & |
---|
| 187 | / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & |
---|
| 188 | - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) |
---|
| 189 | |
---|
| 190 | ELSE |
---|
| 191 | vt = 0 |
---|
| 192 | ENDIF |
---|
| 193 | d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc) |
---|
| 194 | d2=cos(latc)*sin(lon-lonc) |
---|
| 195 | d=MAX(1e-25,sqrt(d1**2+d2**2)) |
---|
| 196 | ulon(ij+u_lup,l)=vt*d1/d |
---|
| 197 | ulat(ij+u_lup,l)=vt*d2/d |
---|
| 198 | |
---|
| 199 | |
---|
| 200 | |
---|
| 201 | |
---|
| 202 | CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) |
---|
| 203 | zz=0.5*(z(ij,l)+z(ij+t_ldown,l)) |
---|
| 204 | r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) |
---|
| 205 | IF (zz<=zt) THEN |
---|
| 206 | vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & |
---|
| 207 | / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & |
---|
| 208 | - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) |
---|
| 209 | |
---|
| 210 | ELSE |
---|
| 211 | vt = 0 |
---|
| 212 | ENDIF |
---|
| 213 | d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc) |
---|
| 214 | d2=cos(latc)*sin(lon-lonc) |
---|
| 215 | d=MAX(1e-25,sqrt(d1**2+d2**2)) |
---|
| 216 | ulon(ij+u_ldown,l)=vt*d1/d |
---|
| 217 | ulat(ij+u_ldown,l)=vt*d2/d |
---|
| 218 | |
---|
| 219 | |
---|
| 220 | |
---|
| 221 | CALL xyz2lonlat(xyz_i(ij,:),lon,lat) |
---|
| 222 | zz=z(ij,l) |
---|
| 223 | r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) |
---|
| 224 | IF (zz<=zt) THEN |
---|
| 225 | vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & |
---|
| 226 | / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & |
---|
| 227 | - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) |
---|
| 228 | ELSE |
---|
| 229 | vt = 0 |
---|
| 230 | ENDIF |
---|
| 231 | d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc) |
---|
| 232 | d2=cos(latc)*sin(lon-lonc) |
---|
| 233 | d=MAX(1e-25,sqrt(d1**2+d2**2)) |
---|
| 234 | ! ulon(ij+u_ldown,l)=vt*d1/d |
---|
| 235 | ! ulat(ij+u_ldown,l)=vt*d2/d |
---|
[114] | 236 | ! out_i(ij,l)=sqrt((vt*d1/d)**2+(vt*d2/d)**2) |
---|
[113] | 237 | ENDDO |
---|
| 238 | ENDDO |
---|
| 239 | ENDDO |
---|
| 240 | |
---|
| 241 | CALL compute_wind_perp_from_lonlat_compound(ulon, ulat, ue) |
---|
[114] | 242 | CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) |
---|
| 243 | |
---|
[113] | 244 | END SUBROUTINE compute_etat0_dcmip5 |
---|
| 245 | |
---|
| 246 | SUBROUTINE compute_z(pmodel,ps,r,z) |
---|
| 247 | USE icosa |
---|
| 248 | IMPLICIT NONE |
---|
| 249 | REAL(rstd),PARAMETER :: epsilon0=2e-13 |
---|
| 250 | REAL(rstd),INTENT(IN) :: pmodel |
---|
| 251 | REAL(rstd),INTENT(IN) :: ps |
---|
| 252 | REAL(rstd),INTENT(IN) :: r |
---|
| 253 | REAL(rstd),INTENT(OUT) :: z |
---|
| 254 | |
---|
| 255 | REAL(rstd) :: Tv0,Tvt,pt |
---|
| 256 | REAL(rstd) :: pave, pp, dpdz |
---|
| 257 | REAL(rstd) :: znp1 |
---|
| 258 | REAL(rstd) :: epsilon |
---|
| 259 | REAL(rstd) :: p |
---|
| 260 | INTEGER :: n |
---|
| 261 | |
---|
| 262 | Tv0=T0/(1+0.608*q0) |
---|
| 263 | Tvt=Tv0-Gamma*zt |
---|
| 264 | pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma)) |
---|
| 265 | |
---|
| 266 | IF (pmodel>pt) THEN |
---|
| 267 | z=Tv0/Gamma*(1-(pmodel/ps)**(Rd*Gamma/g)) |
---|
| 268 | ELSE |
---|
| 269 | z=zt+Rd*Tvt/g*log(Pt/pmodel) |
---|
| 270 | ENDIF |
---|
| 271 | |
---|
| 272 | ! PRINT*,pmodel,pt,z,zt |
---|
| 273 | IF (z>zt .OR. r>1000000) return |
---|
| 274 | |
---|
| 275 | epsilon=1 |
---|
| 276 | n=0 |
---|
| 277 | |
---|
| 278 | DO WHILE(epsilon>epsilon0 .AND. n<20) |
---|
| 279 | |
---|
| 280 | IF (z<zt) THEN |
---|
| 281 | pave=pb*((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma)) |
---|
| 282 | pp= -Deltap * exp(-(r/rp)**1.5-(z/zp)**2) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)) |
---|
| 283 | ELSE |
---|
| 284 | pave=pt*exp(g*(zt-z)/(Rd*Tvt)) |
---|
| 285 | pp=0 |
---|
| 286 | ENDIF |
---|
| 287 | p=pave+pp |
---|
| 288 | |
---|
| 289 | dpdz = 2*Deltap*z/zp**2 * exp(-(r/rp)**1.5) * exp(-(z/zp)**2) * ((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma)) & |
---|
| 290 | - g/(Rd*Tv0) * (preff - Deltap*exp(-(r/rp)**1.5) * exp(-(z/zp)**2)) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)-1) |
---|
| 291 | |
---|
| 292 | znp1 = z - ( pmodel - p ) / (-dpdz) |
---|
| 293 | |
---|
| 294 | epsilon=ABS( (znp1-z)/znp1 ) |
---|
| 295 | ! PRINT *,"epsilon",epsilon,r,z,znp1 |
---|
| 296 | ! PRINT *,"----",n,pmodel,ps,pave,dpdz |
---|
| 297 | z=znp1 |
---|
| 298 | n=n+1 |
---|
| 299 | ENDDO |
---|
| 300 | |
---|
| 301 | END SUBROUTINE compute_z |
---|
| 302 | |
---|
| 303 | END MODULE etat0_dcmip5_mod |
---|