![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 203 for codes/icosagcm/trunk/src/etat0_jablonowsky06.f90
- Timestamp:
- 07/09/14 00:58:30 (10 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/etat0_jablonowsky06.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/etat0_jablonowsky06.f90
r201 r203 11 11 REAL(rstd),PARAMETER :: Gamma=0.005 12 12 REAL(rstd),PARAMETER :: up0=1 13 PUBLIC test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06, compute_etat0_new 13 14 PUBLIC compute_etat0 15 14 16 CONTAINS 15 17 16 SUBROUTINE test_etat0_jablonowsky06 17 USE icosa 18 USE kinetic_mod 19 USE pression_mod 20 USE exner_mod 21 USE geopotential_mod 22 USE vorticity_mod 23 IMPLICIT NONE 24 TYPE(t_field),POINTER,SAVE :: f_ps(:) 25 TYPE(t_field),POINTER,SAVE :: f_phis(:) 26 TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:) 27 TYPE(t_field),POINTER,SAVE :: f_u(:) 28 TYPE(t_field),POINTER,SAVE :: f_Ki(:) 29 TYPE(t_field),POINTER,SAVE :: f_temp(:) 30 TYPE(t_field),POINTER,SAVE :: f_p(:) 31 TYPE(t_field),POINTER,SAVE :: f_pks(:) 32 TYPE(t_field),POINTER,SAVE :: f_pk(:) 33 TYPE(t_field),POINTER,SAVE :: f_phi(:) 34 TYPE(t_field),POINTER,SAVE :: f_vort(:) 35 TYPE(t_field),POINTER,SAVE :: f_q(:) 36 37 REAL(rstd),POINTER :: Ki(:,:) 38 REAL(rstd),POINTER :: temp(:) 39 INTEGER :: ind 40 41 42 CALL allocate_field(f_ps,field_t,type_real) 43 CALL allocate_field(f_phis,field_t,type_real) 44 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 45 CALL allocate_field(f_u,field_u,type_real,llm) 46 CALL allocate_field(f_p,field_t,type_real,llm+1) 47 CALL allocate_field(f_Ki,field_t,type_real,llm) 48 CALL allocate_field(f_pks,field_t,type_real) 49 CALL allocate_field(f_pk,field_t,type_real,llm) 50 CALL allocate_field(f_phi,field_t,type_real,llm) 51 CALL allocate_field(f_temp,field_t,type_real) 52 CALL allocate_field(f_vort,field_z,type_real,llm) 53 54 CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 55 56 CALL kinetic(f_u,f_Ki) 57 CALL vorticity(f_u,f_vort) 58 59 CALL pression(f_ps,f_p) 60 CALL exner(f_ps,f_p,f_pks,f_pk) 61 CALL geopotential(f_phis,f_pks,f_pk,f_theta_rhodz,f_phi) 62 63 CALL writefield('ps',f_ps) 64 CALL writefield('phis',f_phis) 65 CALL writefield('theta',f_theta_rhodz) 66 CALL writefield('f_phi',f_phi) 67 68 CALL writefield('Ki',f_Ki) 69 CALL writefield('vort',f_vort) 70 71 END SUBROUTINE test_etat0_jablonowsky06 72 73 74 SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 75 USE icosa 76 IMPLICIT NONE 77 TYPE(t_field),POINTER :: f_ps(:) 78 TYPE(t_field),POINTER :: f_phis(:) 79 TYPE(t_field),POINTER :: f_theta_rhodz(:) 80 TYPE(t_field),POINTER :: f_u(:) 81 TYPE(t_field),POINTER :: f_q(:) 82 83 REAL(rstd),POINTER :: ps(:) 84 REAL(rstd),POINTER :: phis(:) 85 REAL(rstd),POINTER :: theta_rhodz(:,:) 86 REAL(rstd),POINTER :: u(:,:) 87 REAL(rstd),POINTER :: q(:,:,:) 88 INTEGER :: ind 89 90 DO ind=1,ndomain 91 IF (.NOT. assigned_domain(ind)) CYCLE 92 CALL swap_dimensions(ind) 93 CALL swap_geometry(ind) 94 ps=f_ps(ind) 95 phis=f_phis(ind) 96 theta_rhodz=f_theta_rhodz(ind) 97 u=f_u(ind) 98 q=f_q(ind) 99 q=0 100 CALL compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u) 101 ENDDO 102 103 END SUBROUTINE etat0 104 105 SUBROUTINE compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u) 106 USE icosa 107 USE disvert_mod 108 USE pression_mod 109 USE exner_mod 110 USE geopotential_mod 111 USE theta2theta_rhodz_mod 112 IMPLICIT NONE 113 REAL(rstd),INTENT(OUT) :: ps(iim*jjm) 114 REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 115 REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 116 REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 117 118 INTEGER :: i,j,l,ij 119 REAL(rstd) :: theta(iim*jjm,llm) 120 REAL(rstd) :: eta(llm) 121 REAL(rstd) :: etav(llm) 122 REAL(rstd) :: etas, etavs 123 REAL(rstd) :: lon,lat 124 REAL(rstd) :: ulon(3) 125 REAL(rstd) :: ep(3), norm_ep 126 REAL(rstd) :: Tave, T 127 REAL(rstd) :: phis_ave 128 REAL(rstd) :: V0(3) 129 REAL(rstd) :: r2 130 REAL(rstd) :: utot 131 REAL(rstd) :: lonx,latx 132 133 DO l=1,llm 134 eta(l)= 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) ) 135 etav(l)=(eta(l)-eta0)*Pi/2 136 ENDDO 137 etas=ap(1)+bp(1) 138 etavs=(etas-eta0)*Pi/2 139 140 DO j=jj_begin,jj_end 141 DO i=ii_begin,ii_end 142 ij=(j-1)*iim+i 143 ps(ij)=ps0 144 ENDDO 145 ENDDO 146 147 148 CALL lonlat2xyz(Pi/9,2*Pi/9,V0) 149 150 u(:,:)=1e10 151 DO l=1,llm 152 DO j=jj_begin,jj_end 153 DO i=ii_begin,ii_end 154 ij=(j-1)*iim+i 155 156 CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) 157 CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep) 158 r2=(asin(sqrt(sum(ep*ep))))**2 159 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 160 161 ulon(1) = -sin(lon) * utot 162 ulon(2) = cos(lon) * utot 163 ulon(3) = 0 * utot 164 165 166 CALL cross_product2(xyz_v(ij+z_rdown,:)/radius,xyz_v(ij+z_rup,:)/radius,ep) 167 norm_ep=sqrt(sum(ep(:)**2)) 168 IF (norm_ep>1e-30) THEN 169 ep=-ep/norm_ep*ne(ij,right) 170 u(ij+u_right,l)=sum(ep(:)*ulon(:)) 171 ENDIF 172 173 174 175 CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) 176 CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep) 177 r2=(asin(sqrt(sum(ep*ep))))**2 178 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 179 ulon(1) = -sin(lon) * utot 180 ulon(2) = cos(lon) * utot 181 ulon(3) = 0 * utot 182 183 184 CALL cross_product2(xyz_v(ij+z_up,:)/radius,xyz_v(ij+z_lup,:)/radius,ep) 185 norm_ep=sqrt(sum(ep(:)**2)) 186 IF (norm_ep>1e-30) THEN 187 ep=-ep/norm_ep*ne(ij,lup) 188 u(ij+u_lup,l)=sum(ep(:)*ulon(:)) 189 ENDIF 190 191 192 193 194 195 196 CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) 197 CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep) 198 r2=(asin(sqrt(sum(ep*ep))))**2 199 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 200 ulon(1) = -sin(lon) * utot 201 ulon(2) = cos(lon) * utot 202 ulon(3) = 0 * utot 203 204 205 CALL cross_product2(xyz_v(ij+z_ldown,:)/radius,xyz_v(ij+z_down,:)/radius,ep) 206 norm_ep=sqrt(sum(ep(:)**2)) 207 IF (norm_ep>1e-30) THEN 208 ep=-ep/norm_ep*ne(ij,ldown) 209 u(ij+u_ldown,l)=sum(ep(:)*ulon(:)) 210 ENDIF 211 212 ENDDO 213 ENDDO 214 ENDDO 215 216 217 DO l=1,llm 218 Tave=T0*eta(l)**(Rd*Gamma/g) 219 IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 220 DO j=jj_begin,jj_end 221 DO i=ii_begin,ii_end 222 ij=(j-1)*iim+i 223 CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 224 225 T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 & 226 * ( (-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 & 227 + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega) 228 229 theta(ij,l)=T*eta(l)**(-kappa) 230 231 ENDDO 232 ENDDO 233 ENDDO 234 235 236 phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) 237 DO j=jj_begin,jj_end 238 DO i=ii_begin,ii_end 239 ij=(j-1)*iim+i 240 CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 241 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 & 242 +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega ) 243 ENDDO 244 ENDDO 245 246 CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) 247 248 END SUBROUTINE compute_etat0_jablonowsky06 249 250 SUBROUTINE compute_etat0_new(ngrid,lat,lon, phis, ps, temp, ulon, ulat) 18 SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, temp, ulon, ulat) 251 19 USE disvert_mod 252 20 IMPLICIT NONE … … 283 51 IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 284 52 DO ij=1,ngrid 285 r2 = arc(Pi/9 ,2*Pi/9, lon(ij),lat(ij))**253 r2 = arc(Pi/9.,2.*Pi/9., lon(ij),lat(ij))**2 286 54 temp(ij,l) = Tave + 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 & 287 55 * ( (-2*sin(lat(ij))**6*(cos(lat(ij))**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 & … … 293 61 294 62 295 END SUBROUTINE compute_etat0 _new63 END SUBROUTINE compute_etat0 296 64 297 65 END MODULE etat0_jablonowsky06_mod
Note: See TracChangeset
for help on using the changeset viewer.
