| 1 | MODULE etat0_dcmip41_mod |
|---|
| 2 | USE icosa |
|---|
| 3 | PRIVATE |
|---|
| 4 | REAL(rstd),PARAMETER :: eta0=0.252 |
|---|
| 5 | REAL(rstd),PARAMETER :: etat=0.2 |
|---|
| 6 | REAL(rstd),PARAMETER :: ps0=1e5 |
|---|
| 7 | REAL(rstd),PARAMETER :: u0=35 |
|---|
| 8 | REAL(rstd),PARAMETER :: T0=288 |
|---|
| 9 | REAL(rstd),PARAMETER :: DeltaT=4.8e5 |
|---|
| 10 | REAL(rstd),PARAMETER :: Rd=287 |
|---|
| 11 | REAL(rstd),PARAMETER :: Gamma=0.005 |
|---|
| 12 | REAL(rstd),PARAMETER :: up0=1 |
|---|
| 13 | REAL(rstd) :: lonc |
|---|
| 14 | REAL(rstd) :: latc |
|---|
| 15 | PUBLIC etat0 |
|---|
| 16 | CONTAINS |
|---|
| 17 | |
|---|
| 18 | |
|---|
| 19 | |
|---|
| 20 | SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) |
|---|
| 21 | USE icosa |
|---|
| 22 | IMPLICIT NONE |
|---|
| 23 | TYPE(t_field),POINTER :: f_ps(:) |
|---|
| 24 | TYPE(t_field),POINTER :: f_phis(:) |
|---|
| 25 | TYPE(t_field),POINTER :: f_theta_rhodz(:) |
|---|
| 26 | TYPE(t_field),POINTER :: f_u(:) |
|---|
| 27 | TYPE(t_field),POINTER :: f_q(:) |
|---|
| 28 | |
|---|
| 29 | REAL(rstd),POINTER :: ps(:) |
|---|
| 30 | REAL(rstd),POINTER :: phis(:) |
|---|
| 31 | REAL(rstd),POINTER :: theta_rhodz(:,:) |
|---|
| 32 | REAL(rstd),POINTER :: u(:,:) |
|---|
| 33 | REAL(rstd),POINTER :: q(:,:,:) |
|---|
| 34 | INTEGER :: ind |
|---|
| 35 | |
|---|
| 36 | DO ind=1,ndomain |
|---|
| 37 | CALL swap_dimensions(ind) |
|---|
| 38 | CALL swap_geometry(ind) |
|---|
| 39 | ps=f_ps(ind) |
|---|
| 40 | phis=f_phis(ind) |
|---|
| 41 | theta_rhodz=f_theta_rhodz(ind) |
|---|
| 42 | u=f_u(ind) |
|---|
| 43 | q=f_q(ind) |
|---|
| 44 | CALL compute_etat0_dcmip41(ps, phis, theta_rhodz, u, q) |
|---|
| 45 | ENDDO |
|---|
| 46 | |
|---|
| 47 | END SUBROUTINE etat0 |
|---|
| 48 | |
|---|
| 49 | SUBROUTINE compute_etat0_dcmip41(ps, phis, theta_rhodz, u, q) |
|---|
| 50 | USE icosa |
|---|
| 51 | USE disvert_mod |
|---|
| 52 | USE pression_mod |
|---|
| 53 | USE exner_mod |
|---|
| 54 | USE geopotential_mod |
|---|
| 55 | USE theta2theta_rhodz_mod |
|---|
| 56 | IMPLICIT NONE |
|---|
| 57 | REAL(rstd),INTENT(OUT) :: ps(iim*jjm) |
|---|
| 58 | REAL(rstd),INTENT(OUT) :: phis(iim*jjm) |
|---|
| 59 | REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) |
|---|
| 60 | REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) |
|---|
| 61 | REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) |
|---|
| 62 | |
|---|
| 63 | INTEGER :: i,j,l,ij |
|---|
| 64 | REAL(rstd) :: theta(iim*jjm,llm) |
|---|
| 65 | REAL(rstd) :: Y(iim*jjm,llm) |
|---|
| 66 | REAL(rstd) :: vort |
|---|
| 67 | REAL(rstd) :: eta(llm) |
|---|
| 68 | REAL(rstd) :: etav(llm) |
|---|
| 69 | REAL(rstd) :: etas, etavs |
|---|
| 70 | REAL(rstd) :: lon,lat |
|---|
| 71 | REAL(rstd) :: ulon(3) |
|---|
| 72 | REAL(rstd) :: ep(3), norm_ep |
|---|
| 73 | REAL(rstd) :: Tave, T |
|---|
| 74 | REAL(rstd) :: phis_ave |
|---|
| 75 | REAL(rstd) :: V0(3) |
|---|
| 76 | REAL(rstd) :: r2 |
|---|
| 77 | REAL(rstd) :: utot |
|---|
| 78 | REAL(rstd) :: lonx,latx |
|---|
| 79 | REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r |
|---|
| 80 | lonc=Pi/9 |
|---|
| 81 | latc=2*Pi/9 |
|---|
| 82 | |
|---|
| 83 | DO l=1,llm |
|---|
| 84 | eta(l)= 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) ) |
|---|
| 85 | etav(l)=(eta(l)-eta0)*Pi/2 |
|---|
| 86 | ENDDO |
|---|
| 87 | etas=ap(1)/preff+bp(1) |
|---|
| 88 | etavs=(etas-eta0)*Pi/2 |
|---|
| 89 | |
|---|
| 90 | DO j=jj_begin,jj_end |
|---|
| 91 | DO i=ii_begin,ii_end |
|---|
| 92 | ij=(j-1)*iim+i |
|---|
| 93 | ps(ij)=ps0 |
|---|
| 94 | ENDDO |
|---|
| 95 | ENDDO |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | CALL lonlat2xyz(lonc,latc,V0) |
|---|
| 99 | |
|---|
| 100 | u(:,:)=1e10 |
|---|
| 101 | DO l=1,llm |
|---|
| 102 | DO j=jj_begin-1,jj_end+1 |
|---|
| 103 | DO i=ii_begin-1,ii_end+1 |
|---|
| 104 | ij=(j-1)*iim+i |
|---|
| 105 | |
|---|
| 106 | CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) |
|---|
| 107 | K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) |
|---|
| 108 | r=radius*acos(K) |
|---|
| 109 | utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) |
|---|
| 110 | u(ij+u_right,l) = utot * sum(elon_e(ij+u_right,:) * ep_e(ij+u_right,:)) |
|---|
| 111 | |
|---|
| 112 | |
|---|
| 113 | CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) |
|---|
| 114 | K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) |
|---|
| 115 | r=radius*acos(K) |
|---|
| 116 | utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) |
|---|
| 117 | u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:)) |
|---|
| 118 | |
|---|
| 119 | CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) |
|---|
| 120 | K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) |
|---|
| 121 | r=radius*acos(K) |
|---|
| 122 | utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) |
|---|
| 123 | u(ij+u_ldown,l) = utot * sum(elon_e(ij+u_ldown,:) * ep_e(ij+u_ldown,:)) |
|---|
| 124 | |
|---|
| 125 | ENDDO |
|---|
| 126 | ENDDO |
|---|
| 127 | ENDDO |
|---|
| 128 | |
|---|
| 129 | |
|---|
| 130 | DO l=1,llm |
|---|
| 131 | Tave=T0*eta(l)**(Rd*Gamma/g) |
|---|
| 132 | IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 |
|---|
| 133 | DO j=jj_begin,jj_end |
|---|
| 134 | DO i=ii_begin,ii_end |
|---|
| 135 | ij=(j-1)*iim+i |
|---|
| 136 | CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) |
|---|
| 137 | |
|---|
| 138 | Y(ij,l)=((-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 & |
|---|
| 139 | + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega) |
|---|
| 140 | T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) |
|---|
| 141 | |
|---|
| 142 | theta(ij,l)=T*eta(l)**(-kappa) |
|---|
| 143 | |
|---|
| 144 | ENDDO |
|---|
| 145 | ENDDO |
|---|
| 146 | ENDDO |
|---|
| 147 | |
|---|
| 148 | |
|---|
| 149 | phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) |
|---|
| 150 | DO j=jj_begin,jj_end |
|---|
| 151 | DO i=ii_begin,ii_end |
|---|
| 152 | ij=(j-1)*iim+i |
|---|
| 153 | CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) |
|---|
| 154 | phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5 & |
|---|
| 155 | +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega ) |
|---|
| 156 | ! phis(ij)=phis_ave+u0*cos(etavs)**1.5 |
|---|
| 157 | |
|---|
| 158 | ENDDO |
|---|
| 159 | ENDDO |
|---|
| 160 | |
|---|
| 161 | CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) |
|---|
| 162 | |
|---|
| 163 | q(:,:,1)=theta(:,:) |
|---|
| 164 | |
|---|
| 165 | |
|---|
| 166 | DO l=1,llm |
|---|
| 167 | dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1) |
|---|
| 168 | IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa) & |
|---|
| 169 | + kappa * (etat-eta(l))**5 * eta(l)**(-kappa-1)) |
|---|
| 170 | DO j=jj_begin,jj_end |
|---|
| 171 | DO i=ii_begin,ii_end |
|---|
| 172 | ij=(j-1)*iim+i |
|---|
| 173 | CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) |
|---|
| 174 | dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) & |
|---|
| 175 | + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l) & |
|---|
| 176 | - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)& |
|---|
| 177 | - 9./8. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l)) & |
|---|
| 178 | * (-2*sin(lat)**6*(cos(lat)**2+1./3.)+10./63.) |
|---|
| 179 | dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5 & |
|---|
| 180 | *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) & |
|---|
| 181 | + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat))) |
|---|
| 182 | |
|---|
| 183 | duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l)) |
|---|
| 184 | |
|---|
| 185 | K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) |
|---|
| 186 | r=radius*acos(K) |
|---|
| 187 | vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2) & |
|---|
| 188 | + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) & |
|---|
| 189 | -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2))) |
|---|
| 190 | q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta)) |
|---|
| 191 | ENDDO |
|---|
| 192 | ENDDO |
|---|
| 193 | ENDDO |
|---|
| 194 | |
|---|
| 195 | |
|---|
| 196 | END SUBROUTINE compute_etat0_dcmip41 |
|---|
| 197 | |
|---|
| 198 | END MODULE etat0_dcmip41_mod |
|---|
![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
