Changeset 741 for codes/icosagcm/devel
- Timestamp:
- 09/28/18 12:59:46 (6 years ago)
- Location:
- codes/icosagcm/devel/src
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/initial/etat0_venus.f90
r531 r741 3 3 IMPLICIT NONE 4 4 PRIVATE 5 SAVE 5 6 6 7 TYPE(t_field),POINTER :: f_temp_eq( :) 7 8 TYPE(t_field),POINTER :: f_temp(:) ! buffer used for physics 8 9 REAL(rstd), SAVE :: kfrict 10 !$OMP THREADPRIVATE(kfrict) 11 12 PUBLIC :: etat0, init_physics, physics 9 REAL(rstd), ALLOCATABLE :: temp_eq_packed(:,:) 10 11 REAL(rstd) :: kfrict, kv 12 !$OMP THREADPRIVATE(kfrict, kv) 13 14 REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012 15 16 PUBLIC :: etat0, init_physics, full_physics 13 17 14 18 CONTAINS … … 16 20 !-------------------------------- "Physics" ---------------------------------------- 17 21 18 SUBROUTINE init_physics 22 SUBROUTINE init_physics_old 19 23 USE getin_mod 20 IMPLICIT NONE21 24 REAL(rstd),POINTER :: temp(:,:) 22 25 REAL(rstd) :: friction_time 23 26 INTEGER :: ind 24 27 28 kv=0.15 ! vertical turbulent viscosity 29 CALL getin('venus_diffusion',kv) 25 30 friction_time=86400. !friction 26 CALL getin(' friction_time',friction_time)31 CALL getin('venus_friction_time',friction_time) 27 32 kfrict=1./friction_time 28 33 29 CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use 34 CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use by physics 30 35 31 36 PRINT *, 'Initializing Temp_eq (venus)' … … 36 41 CALL swap_geometry(ind) 37 42 temp=f_temp_eq(ind) 38 CALL compute_temp_ref( temp, .TRUE.) ! FIXMEWith meridional gradient39 ENDDO 40 END SUBROUTINE init_physics 43 CALL compute_temp_ref(.TRUE.,iim*jjm,lat_i, temp) ! With meridional gradient 44 ENDDO 45 END SUBROUTINE init_physics_old 41 46 42 47 SUBROUTINE physics(f_ps,f_theta_rhodz,f_u) 43 USE icosa44 48 USE theta2theta_rhodz_mod 45 IMPLICIT NONE46 49 TYPE(t_field),POINTER :: f_theta_rhodz(:) 47 50 TYPE(t_field),POINTER :: f_u(:) … … 66 69 67 70 SUBROUTINE compute_physics(temp_eq, temp, u) 68 USE icosa69 USE theta2theta_rhodz_mod70 IMPLICIT NONE71 71 REAL(rstd),INTENT(IN) :: temp_eq(iim*jjm,llm) 72 72 REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm) 73 73 REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm) 74 REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 201275 74 INTEGER :: i,j,l,ij 76 75 DO l=1,llm … … 85 84 END SUBROUTINE compute_physics 86 85 86 !----------------------- Re-implementation using physics_interface_mod ----------------- 87 88 SUBROUTINE init_physics 89 USE getin_mod 90 USE physics_interface_mod 91 REAL(rstd),POINTER :: temp(:,:) 92 REAL(rstd) :: friction_time 93 INTEGER :: ngrid 94 95 kv=0.15 ! vertical turbulent viscosity 96 CALL getin('venus_diffusion',kv) 97 friction_time=86400. !friction 98 CALL getin('venus_friction_time',friction_time) 99 kfrict=1./friction_time 100 101 ngrid = physics_inout%ngrid 102 ALLOCATE(temp_eq_packed(ngrid,llm)) 103 CALL compute_temp_ref(.TRUE.,ngrid,physics_inout%lat, temp_eq_packed) ! With meridional gradient 104 END SUBROUTINE init_physics 105 106 SUBROUTINE full_physics 107 USE physics_interface_mod 108 CALL compute_physics_column(physics_inout%ngrid, physics_inout%dt_phys, & 109 physics_inout%p, physics_inout%geopot, & 110 physics_inout%Temp, physics_inout%ulon, physics_inout%ulat, & 111 physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat) 112 END SUBROUTINE full_physics 113 114 SUBROUTINE compute_physics_column(ngrid,dt_phys,p,geopot,Temp,u,v,dTemp,du,dv) 115 USE earth_const, only : g 116 INTEGER, INTENT(IN) :: ngrid 117 REAL(rstd),INTENT(IN) :: dt_phys, & 118 p(ngrid,llm+1), geopot(ngrid,llm+1), & 119 Temp(ngrid,llm), u(ngrid,llm), v(ngrid,llm) 120 REAL(rstd),INTENT(OUT) :: dTemp(ngrid,llm), du(ngrid,llm), dv(ngrid,llm) 121 REAL(rstd) :: mass(ngrid,llm), & ! rho.dz 122 A(ngrid,llm+1), & ! off-diagonal coefficients 123 B(ngrid,llm), & ! diagonal coefficients 124 Ru(ngrid,llm), Rv(ngrid,llm), & ! right-hand-sides 125 C(ngrid,llm), & ! LU factorization (Thomas algorithm) 126 xu(ngrid,llm), xv(ngrid, llm) ! solution 127 REAL(rstd) :: rho, X_ij 128 INTEGER :: l,ij 129 ! Vertical diffusion : rho.du/dt = d/dz (rho*kappa*du/dz) 130 ! rho.dz = mass.deta 131 ! => mass.du/dt = d/deta ((rho^2 kappa/m)du/deta) 132 ! Backward Euler : (M-S)u_new = M.u_old 133 ! with M.u = mass(l)*u(l) 134 ! and S.u = A(l+1)*(u(l+1)-u(l)) - A(l)*(u(l)-u(l-1)) 135 ! => solve -A(l)u(l-1) + B(l)u(l) - A(l+1)u(l+1) = Ru(l) using Thomas algorithm 136 ! with A=tau*kappa*rho^2/m and B(l) = mass(l)+A(l)+A(l+1) 137 DO l=1,llm 138 !DIR$ SIMD 139 DO ij=1,ngrid 140 mass(ij,l) = (p(ij,l)-p(ij,l+1))*(1./g) 141 rho = (p(ij,l)-p(ij,l+1))/(geopot(ij,l+1)-geopot(ij,l)) 142 ! A = kappa.tau.rho^2/m 143 A(ij,l) = (kv*dt_phys)*(rho**2)/mass(ij,l) 144 Ru(ij,l) = mass(ij,l)*u(ij,l) 145 Rv(ij,l) = mass(ij,l)*v(ij,l) 146 ENDDO 147 ENDDO 148 A(:,llm+1)=0. 149 DO l=llm,2,-1 150 !DIR$ SIMD 151 DO ij=1,ngrid 152 A(ij,l) = .5*(A(ij,l)+A(ij,l-1)) ! average A to interfaces 153 ENDDO 154 ENDDO 155 A(:,1)=0. 156 DO l=1,llm 157 !DIR$ SIMD 158 DO ij=1,ngrid 159 B(ij,l) = mass(ij,l)+A(ij,l)+A(ij,l+1) 160 ENDDO 161 ENDDO 162 163 ! Solve -A(l)x(l-1) + B(l)x(l) - A(l+1)x(l+1) = R(l) using Thomas algorithm 164 ! Forward sweep : 165 ! C(0)=0, C(l) = -A(l+1) / (B(l)+A(l)C(l-1)), 166 ! D(0)=0, D(l) = (R(l)+A(l)D(l-1)) / (B(l)+A(l)C(l-1)) 167 !DIR$ SIMD 168 DO ij=1,ngrid 169 X_ij = 1./B(ij,1) 170 C(ij,2) = -A(ij,2) * X_ij 171 xu(ij,1) = Ru(ij,1) * X_ij 172 xv(ij,1) = Rv(ij,1) * X_ij 173 END DO 174 DO l = 2,llm 175 !DIR$ SIMD 176 DO ij=1,ngrid 177 X_ij = 1./( B(ij,l) + A(ij,l)*C(ij,l) ) 178 C(ij,l+1) = -A(ij,l+1) * X_ij ! zero for l=llm 179 xu(ij,l) = (Ru(ij,l)+A(ij,l)*xu(ij,l-1)) * X_ij 180 xv(ij,l) = (Rv(ij,l)+A(ij,l)*xv(ij,l-1)) * X_ij 181 END DO 182 END DO 183 ! Back substitution : 184 ! x(i) = D(i)-C(i+1)x(i+1), x(llm)=D(llm) 185 !DIR$ SIMD 186 DO ij=1,ngrid 187 ! top layer l=llm 188 du(ij,llm) = (xu(ij,llm)-u(ij,llm))/dt_phys 189 dv(ij,llm) = (xv(ij,llm)-v(ij,llm))/dt_phys 190 END DO 191 ! Back substitution at lower layers 192 DO l = llm-1,1,-1 193 !DIR$ SIMD 194 DO ij=1,ngrid 195 xu(ij,l) = xu(ij,l) - C(ij,l+1)*xu(ij,l+1) 196 xv(ij,l) = xv(ij,l) - C(ij,l+1)*xv(ij,l+1) 197 du(ij,l) = (xu(ij,l)-u(ij,l))/dt_phys 198 dv(ij,l) = (xv(ij,l)-v(ij,l))/dt_phys 199 END DO 200 END DO 201 ! bottom friction + thermal relaxation 202 du(:,1)=du(:,1)-kfrict*u(:,1) 203 dTemp(:,:) = (1./tauCLee)*(temp_eq_packed(:,:)-temp(:,:)) 204 END SUBROUTINE compute_physics_column 205 87 206 !----------------------------- Initialize to T_eq -------------------------------------- 88 207 89 208 SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 90 USE icosa91 209 USE theta2theta_rhodz_mod 92 IMPLICIT NONE93 210 TYPE(t_field),POINTER :: f_ps(:) 94 211 TYPE(t_field),POINTER :: f_phis(:) … … 103 220 REAL(rstd),POINTER :: u(:,:) 104 221 REAL(rstd),POINTER :: q(:,:,:) 105 REAL(rstd) :: lat(iim*jjm) ! latitude106 REAL(rstd) :: pplay(iim*jjm, llm) ! pressure at full layers107 222 INTEGER :: ind 108 223 109 CALL allocate_field(f_temp,field_t,type_real,llm )224 CALL allocate_field(f_temp,field_t,type_real,llm, name='temp_buf_venus') 110 225 DO ind=1,ndomain 111 226 IF (.NOT. assigned_domain(ind)) CYCLE … … 121 236 q(:,:,:)=1e2 122 237 temp=f_temp(ind) 123 CALL compute_temp_ref( temp, .FALSE.) ! Without meridional gradient238 CALL compute_temp_ref(.FALSE., iim*jjm, lat_i, temp) ! Without meridional gradient 124 239 ENDDO 125 240 CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) … … 130 245 !------------------------- Compute reference temperature field ------------------------ 131 246 132 SUBROUTINE compute_temp_ref(theta_eq, gradient) 133 USE icosa 247 SUBROUTINE compute_temp_ref(gradient,ngrid,lat, temp_eq) 134 248 USE disvert_mod 135 249 USE omp_para 136 250 USE math_const 137 IMPLICIT NONE138 REAL(rstd), INTENT(OUT) :: theta_eq(iim*jjm,llm)139 251 LOGICAL, INTENT(IN) :: gradient 140 REAL(rstd) :: clat(iim*jjm) ! latitude 141 integer :: level 142 REAL :: pressCLee(31), tempCLee(31), dt_epCLee(31), etaCLee(31) 252 INTEGER, INTENT(IN) :: ngrid 253 REAL(rstd), INTENT(IN) :: lat(ngrid) ! latitude 254 REAL(rstd), INTENT(OUT) :: temp_eq(ngrid,llm) 255 256 REAL(rstd) :: clat(ngrid) 257 INTEGER :: level 258 INTEGER, PARAMETER :: nlev=30 259 REAL :: pressCLee(nlev+1), tempCLee(nlev+1), dt_epCLee(nlev+1), etaCLee(nlev+1) 143 260 144 real(rstd) :: lon,lat, pplay, ztemp,zdt,fact 145 logical, save :: firstcall 146 integer :: i,j,ij, l,ll 261 REAL(rstd) :: pplay, ztemp,zdt,fact 262 INTEGER :: ij, l,ll 147 263 148 264 data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, & … … 154 270 2.265e-6/ 155 271 156 DO j=jj_begin-1,jj_end+1157 DO i=ii_begin-1,ii_end+1158 ij=(j-1)*iim+i159 CALL xyz2lonlat(xyz_i(ij,:),lon,lat)160 clat(ij)=cos(lat)161 ENDDO162 ENDDO163 272 164 273 data tempCLee/ 728.187, 715.129, 697.876, 677.284, 654.078, 628.885, & … … 175 284 2.000/ 176 285 177 pressCLee = etaCLee*9.2e6 178 179 DO j=jj_begin-1,jj_end+1 180 DO i=ii_begin-1,ii_end+1 181 ij=(j-1)*iim+i 182 183 DO l = 1, llm 184 185 pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff 186 ! look for largest level such that pressCLee(level) > pplay(ij,l)) 187 ! => pressClee(level+1) < pplay(ij,l) < pressClee(level) 188 level = 1 189 DO ll = 1, 30 ! 30 data levels 190 IF(pressCLee(ll) > pplay) THEN 191 level = ll 192 END IF 193 END DO 194 195 ! interpolate between level and level+1 196 ! interpolation is linear in log(pressure) 197 fact = ( log10(pplay)-log10(pressCLee(level))) & 198 /( log10(pressCLee(level+1))-log10(pressCLee(level)) ) 199 ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact 200 zdt = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact 201 202 IF(gradient) THEN 203 theta_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.) 204 ELSE 205 theta_eq(ij,l) = ztemp 286 pressCLee(:) = etaCLee(:)*9.2e6 287 clat(:)=COS(lat(:)) 288 289 DO ij=1,ngrid 290 DO l = 1, llm 291 292 pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff 293 ! look for largest level such that pressCLee(level) > pplay(ij,l)) 294 ! => pressClee(level+1) < pplay(ij,l) < pressClee(level) 295 level = 1 296 DO ll = 1, nlev ! nlev data levels 297 IF(pressCLee(ll) > pplay) THEN 298 level = ll 206 299 END IF 207 300 END DO 301 302 ! interpolate between level and level+1 303 ! interpolation is linear in log(pressure) 304 fact = ( log10(pplay)-log10(pressCLee(level))) & 305 /( log10(pressCLee(level+1))-log10(pressCLee(level)) ) 306 ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact 307 zdt = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact 308 309 IF(gradient) THEN 310 temp_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.) 311 ELSE 312 temp_eq(ij,l) = ztemp 313 END IF 208 314 END DO 209 315 END DO -
codes/icosagcm/devel/src/physics/physics.f90
r739 r741 12 12 phys_lmdz_generic=21, phys_external=22 13 13 INTEGER :: phys_type 14 TYPE(t_field),POINTER,SAVE :: f_extra_physics_2D(:), f_extra_physics_3D(:)15 14 TYPE(t_field),POINTER,SAVE :: f_dulon(:), f_dulat(:) 16 15 TYPE(t_field),POINTER,SAVE :: f_ulon(:), f_ulat(:) … … 50 49 CASE ('held_suarez') 51 50 phys_type = phys_HS94 52 CASE ('Lebonnois2012')53 phys_type = phys_LB201254 CALL init_phys_venus55 51 CASE ('phys_lmdz_generic') 56 52 phys_type=phys_lmdz_generic … … 79 75 phys_type = phys_DCMIP2016 80 76 CALL init_physics_dcmip2016 77 CASE ('Lebonnois2012') 78 phys_type = phys_LB2012 79 CALL init_phys_venus 81 80 CASE DEFAULT 82 81 IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',& … … 129 128 END SUBROUTINE add_du_phys 130 129 131 SUBROUTINE physics(it,f_phis, f_ ps, f_theta_rhodz, f_ue, f_wflux, f_q)130 SUBROUTINE physics(it,f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q) 132 131 USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics 133 132 USE physics_external_mod, ONLY : physics_external => physics … … 135 134 USE physics_dcmip2016_mod, ONLY : write_physics_dcmip2016 => write_physics 136 135 USE etat0_heldsz_mod 137 USE etat0_venus_mod, ONLY : phys_venus => physics138 136 INTEGER, INTENT(IN) :: it 139 137 TYPE(t_field),POINTER :: f_phis(:) 138 TYPE(t_field),POINTER :: f_geopot(:) 140 139 TYPE(t_field),POINTER :: f_ps(:) 141 140 TYPE(t_field),POINTER :: f_theta_rhodz(:) … … 161 160 CASE (phys_external) 162 161 CALL physics_external(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q) 163 CASE(phys_LB2012)164 CALL phys_venus(f_ps,f_theta_rhodz,f_ue)165 162 CASE DEFAULT 166 CALL physics_column(it, f_phis, f_ ps, f_theta_rhodz, f_ue, f_q)163 CALL physics_column(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_q) 167 164 END SELECT 168 165 … … 197 194 END SUBROUTINE write_physics_tendencies 198 195 199 SUBROUTINE physics_column(it, f_phis, f_ ps, f_theta_rhodz, f_ue, f_q)196 SUBROUTINE physics_column(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_q) 200 197 USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics 201 198 USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics 199 USE etat0_venus_mod, ONLY : full_physics_venus=>full_physics 202 200 USE theta2theta_rhodz_mod 203 201 USE mpipara 204 202 USE checksum_mod 205 203 TYPE(t_field),POINTER :: f_phis(:) 204 TYPE(t_field),POINTER :: f_geopot(:) 206 205 TYPE(t_field),POINTER :: f_ps(:) 207 206 TYPE(t_field),POINTER :: f_theta_rhodz(:) … … 209 208 TYPE(t_field),POINTER :: f_q(:) 210 209 REAL(rstd),POINTER :: phis(:) 210 REAL(rstd),POINTER :: geopot(:,:) 211 211 REAL(rstd),POINTER :: ps(:) 212 212 REAL(rstd),POINTER :: temp(:,:) … … 228 228 CALL swap_geometry(ind) 229 229 phis=f_phis(ind) 230 geopot=f_geopot(ind) 230 231 ps=f_ps(ind) 231 232 temp=f_temp(ind) … … 236 237 ulon=f_ulon(ind) 237 238 ulat=f_ulat(ind) 238 CALL pack_physics(pack_info(ind), phis, ps, temp, ue, q, p, pk, ulon, ulat)239 CALL pack_physics(pack_info(ind), phis, geopot, ps, temp, ue, q, p, pk, ulon, ulat) 239 240 END DO 240 241 … … 244 245 CASE (phys_DCMIP2016) 245 246 IF (is_omp_level_master) CALL full_physics_dcmip2016 247 CASE(phys_LB2012) 248 IF (is_omp_level_master) CALL full_physics_venus 246 249 CASE DEFAULT 247 250 IF(is_master) PRINT *,'Internal error : illegal value of phys_type', phys_type … … 279 282 END SUBROUTINE physics_column 280 283 281 SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat )284 SUBROUTINE pack_physics(info, phis, geopot, ps, temp, ue, q, p, pk, ulon, ulat ) 282 285 USE wind_mod 283 286 USE pression_mod … … 286 289 TYPE(t_pack_info) :: info 287 290 REAL(rstd) :: phis(iim*jjm) 291 REAL(rstd) :: geopot(iim*jjm,llm+1) 288 292 REAL(rstd) :: ps(iim*jjm) 289 293 REAL(rstd) :: temp(iim*jjm,llm) … … 308 312 IF (is_omp_level_master) THEN 309 313 CALL pack_domain(info, phis, physics_inout%phis) 314 CALL pack_domain(info, geopot, physics_inout%geopot) 310 315 CALL pack_domain(info, p, physics_inout%p) 311 316 CALL pack_domain(info, pk, physics_inout%pk) -
codes/icosagcm/devel/src/physics/physics_interface.f90
r739 r741 11 11 REAL(rstd), DIMENSION(:), POINTER :: Ai, lon, lat, phis 12 12 ! Input, time-dependent 13 REAL(rstd), DIMENSION(:,:), POINTER :: p, pk, Temp, ulon, ulat13 REAL(rstd), DIMENSION(:,:), POINTER :: geopot, p, pk, Temp, ulon, ulat 14 14 REAL(rstd), DIMENSION(:,:,:), POINTER :: q 15 15 ! Output arrays … … 93 93 94 94 ngrid=offset 95 ! Input 95 ! Input 96 96 ALLOCATE(physics_inout%Ai(ngrid)) 97 97 ALLOCATE(physics_inout%lon(ngrid)) … … 99 99 ALLOCATE(physics_inout%phis(ngrid)) 100 100 ALLOCATE(physics_inout%p(ngrid,llm+1)) 101 ALLOCATE(physics_inout%geopot(ngrid,llm+1)) 101 102 ALLOCATE(physics_inout%pk(ngrid,llm)) 102 103 ALLOCATE(physics_inout%Temp(ngrid,llm)) -
codes/icosagcm/devel/src/time/timeloop_gcm.f90
r732 r741 350 350 f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 351 351 CALL enter_profile(id_phys) 352 CALL physics(it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)352 CALL physics(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_u, f_wflux, f_q) 353 353 CALL exit_profile(id_phys) 354 354 CALL check_conserve_detailed(it, AAM_phys, &
Note: See TracChangeset
for help on using the changeset viewer.