Changeset 401 for codes/icosagcm
- Timestamp:
- 06/08/16 01:51:21 (8 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn_gcm.f90
r387 r401 35 35 IF (is_master) PRINT *, 'caldyn_conserv=',def 36 36 37 def='theta' 38 CALL getin('caldyn_thermo',def) 39 SELECT CASE(TRIM(def)) 40 CASE('theta') 41 caldyn_thermo=thermo_theta 42 CASE('entropy') 43 caldyn_thermo=thermo_entropy 44 CASE DEFAULT 45 IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_thermo : <', & 46 TRIM(def),'> options are <theta>, <entropy>' 47 STOP 48 END SELECT 49 37 50 CALL allocate_caldyn 38 51 -
codes/icosagcm/trunk/src/caldyn_kernels_base.f90
r377 r401 37 37 38 38 INTEGER :: i,j,ij,l 39 REAL(rstd) :: p_ik, exner_ik39 REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik 40 40 INTEGER :: ij_omp_begin_ext, ij_omp_end_ext 41 42 41 43 42 CALL trace_start("compute_geopot") … … 47 46 ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 48 47 48 Rd = kappa*cpp 49 49 50 IF(caldyn_eta==eta_mass .AND. .NOT. DEC) THEN 50 51 !!! Compute exner function and geopotential 52 STOP 51 53 DO l = 1,llm 52 54 !DIR$ SIMD … … 89 91 END DO 90 92 ! now pk contains the Lagrange multiplier (pressure) 91 ELSE ! non-Boussinesq, compute geopotential and Exner pressure93 ELSE ! non-Boussinesq, compute pressure, Exner pressure or temperature, then geopotential 92 94 ! uppermost layer 93 95 … … 109 111 END DO 110 112 END IF 111 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 113 112 114 DO l = 1,llm 113 !DIR$ SIMD 114 DO ij=ij_omp_begin_ext,ij_omp_end_ext 115 p_ik = pk(ij,l) 116 exner_ik = cpp * (p_ik/preff) ** kappa 117 pk(ij,l) = exner_ik 118 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 119 ENDDO 115 SELECT CASE(caldyn_thermo) 116 CASE(thermo_theta) 117 !DIR$ SIMD 118 DO ij=ij_omp_begin_ext,ij_omp_end_ext 119 p_ik = pk(ij,l) 120 exner_ik = cpp * (p_ik/preff) ** kappa 121 pk(ij,l) = exner_ik 122 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 123 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 124 ENDDO 125 CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff) 126 !DIR$ SIMD 127 DO ij=ij_omp_begin_ext,ij_omp_end_ext 128 p_ik = pk(ij,l) 129 temp_ik = Treff*exp((theta(ij,l) + Rd*log(p_ik/preff))/cpp) 130 pk(ij,l) = temp_ik 131 ! specific volume v = Rd*T/p = dphi/g/rhodz 132 geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik 133 ENDDO 134 CASE DEFAULT 135 STOP 136 END SELECT 120 137 ENDDO 121 138 END IF -
codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90
r377 r401 401 401 ! from now on pk contains the vertically-averaged geopotential 402 402 pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 403 ENDDO 404 ENDDO 405 403 END DO 404 END DO 406 405 ELSE ! compressible 407 406 408 407 DO l=ll_begin,ll_end 409 !DIR$ SIMD 410 DO ij=ij_begin,ij_end 411 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 412 ENDDO 413 ENDDO 408 SELECT CASE(caldyn_thermo) 409 CASE(thermo_theta) ! vdp = theta.dpi => B = Phi 410 !DIR$ SIMD 411 DO ij=ij_begin,ij_end 412 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 413 END DO 414 CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 415 !DIR$ SIMD 416 DO ij=ij_begin,ij_end 417 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 418 + pk(ij,l)*(cpp-theta(ij,l)) ! pk=temperature, theta=entropy 419 END DO 420 END SELECT 421 END DO 414 422 415 423 END IF ! Boussinesq/compressible -
codes/icosagcm/trunk/src/earth_const.f90
r377 r401 9 9 REAL(rstd),SAVE :: kappa=0.2857143 10 10 REAL(rstd),SAVE :: cpp=1004.70885 11 REAL(rstd),SAVE :: cppv=1004.70885 ! FIXME 12 REAL(rstd),SAVE :: Rv=461.5 13 REAL(rstd),SAVE :: Treff=273. 11 14 REAL(rstd),SAVE :: preff=101325. 12 15 REAL(rstd),SAVE :: pa=50000. … … 15 18 REAL(rstd),SAVE :: gas_constant = 8.3144621 16 19 REAL(rstd),SAVE :: mu ! molar mass of the atmosphere 20 21 INTEGER, PARAMETER,PUBLIC :: thermo_theta=1, thermo_entropy=2, thermo_moist=3 22 INTEGER, PUBLIC :: caldyn_thermo 23 !$OMP THREADPRIVATE(caldyn_thermo) 17 24 18 25 LOGICAL, SAVE :: boussinesq, hydrostatic … … 33 40 CALL getin("cpp",cpp) 34 41 CALL getin("preff",preff) 42 CALL getin("Treff",Treff) 35 43 CALL getin("scale_height",scale_height) 36 44 -
codes/icosagcm/trunk/src/etat0.f90
r388 r401 172 172 CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q) 173 173 ENDIF 174 175 IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1) 176 174 177 ENDDO 175 178 176 IF( TRIM(etat0_type)/='williamson91.6' ) CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)177 178 179 CALL deallocate_field(f_temp) 179 180 180 181 END SUBROUTINE etat0_collocated 182 183 SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset) 184 USE icosa 185 USE pression_mod 186 USE exner_mod 187 USE omp_para 188 IMPLICIT NONE 189 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 190 REAL(rstd),INTENT(IN) :: temp(iim*jjm,llm) 191 REAL(rstd),INTENT(IN) :: q(iim*jjm,llm,nqtot) 192 REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 193 INTEGER,INTENT(IN) :: offset 194 195 REAL(rstd) :: p(iim*jjm,llm+1) 196 REAL(rstd) :: cppd,Rd, mass, p_ij, q_ij,r_ij, chi,nu, entropy, theta 197 INTEGER :: i,j,ij,l 198 199 cppd=cpp 200 Rd=kappa*cppd 201 202 CALL compute_pression(ps,p,offset) 203 ! flush p 204 !$OMP BARRIER 205 DO l = ll_begin, ll_end 206 DO j=jj_begin-offset,jj_end+offset 207 DO i=ii_begin-offset,ii_end+offset 208 ij=(j-1)*iim+i 209 mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass 210 p_ij = .5*(p(ij,l)+p(ij,l+1)) ! pressure at full level 211 SELECT CASE(caldyn_thermo) 212 CASE(thermo_theta) 213 theta = temp(ij,l)*(p_ij/preff)**(-kappa) 214 theta_rhodz(ij,l) = mass * theta 215 CASE(thermo_entropy) 216 nu = log(p_ij/preff) 217 chi = log(temp(ij,l)/Treff) 218 entropy = cppd*chi-Rd*nu 219 theta_rhodz(ij,l) = mass * entropy 220 ! CASE(thermo_moist) 221 ! q_ij=q(ij,l,1) 222 ! r_ij=1.-q_ij 223 ! mass=mass*(1-q_ij) ! dry mass 224 ! nu = log(p_ij/preff) 225 ! chi = log(temp(ij,l)/Treff) 226 ! entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu) 227 ! theta_rhodz(ij,l) = mass * entropy 228 CASE DEFAULT 229 STOP 230 END SELECT 231 ENDDO 232 ENDDO 233 ENDDO 234 !$OMP BARRIER 235 END SUBROUTINE compute_temperature2entropy 181 236 182 237 SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q)
Note: See TracChangeset
for help on using the changeset viewer.