![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 856 for codes/icosagcm/devel/src/initial
- Timestamp:
- 05/06/19 01:49:12 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/initial/etat0_collocated.f90
r848 r856 1 1 MODULE etat0_collocated_mod 2 2 USE icosa 3 USE grid_param 3 4 USE omp_para, ONLY : is_omp_level_master 4 5 USE caldyn_vars_mod, ONLY : hydrostatic … … 6 7 PRIVATE 7 8 8 LOGICAL :: autoinit_mass, autoinit_NH9 9 CHARACTER(len=255),SAVE :: etat0_type 10 !$OMP THREADPRIVATE( autoinit_mass, autoinit_NH,etat0_type)10 !$OMP THREADPRIVATE(etat0_type) 11 11 12 12 PUBLIC :: etat0_type, etat0_collocated … … 36 36 TYPE(t_field),POINTER :: f_q(:) 37 37 38 TYPE(t_field),POINTER,SAVE :: f_temp(:) 38 TYPE(t_field),POINTER,SAVE :: f_temp(:) ! SAVE => all threads see the same f_temp 39 39 REAL(rstd),POINTER :: ps(:) 40 40 REAL(rstd),POINTER :: mass(:,:) … … 64 64 q=f_q(ind) 65 65 66 IF( TRIM(etat0_type)=='williamson91.6' ) THEN 67 CALL compute_etat0_collocated_hex(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q) 68 ELSE 69 CALL compute_etat0_collocated_hex(ps,mass, phis, temp, u, geopot, W, q) 70 ENDIF 71 72 IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1) 73 66 SELECT CASE(grid_type) 67 CASE(grid_ico) 68 CALL compute_etat0_collocated_hex(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q) 69 CASE(grid_unst) 70 CALL compute_etat0_collocated_unst(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q) 71 CASE DEFAULT 72 STOP 'Unexpected value of grid_type encountered in etat0_collocated.' 73 END SELECT 74 74 75 ENDDO 75 76 … … 78 79 END SUBROUTINE etat0_collocated 79 80 80 SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset) 81 USE icosa 82 USE pression_mod 83 USE exner_mod 84 USE omp_para 85 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 86 REAL(rstd),INTENT(IN) :: temp(iim*jjm,llm) 87 REAL(rstd),INTENT(IN) :: q(iim*jjm,llm,nqtot) 88 REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 89 INTEGER,INTENT(IN) :: offset 90 91 REAL(rstd) :: p(iim*jjm,llm+1) 92 REAL(rstd) :: cppd,Rd, mass, p_ij, q_ij,r_ij, chi,nu, entropy, theta 93 INTEGER :: i,j,ij,l 94 95 cppd=cpp 96 Rd=kappa*cppd 97 98 CALL compute_pression(ps,p,offset) 99 ! flush p 100 DO l = ll_begin, ll_end 101 DO j=jj_begin-offset,jj_end+offset 102 DO i=ii_begin-offset,ii_end+offset 103 ij=(j-1)*iim+i 104 mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass 105 p_ij = .5*(p(ij,l)+p(ij,l+1)) ! pressure at full level 106 SELECT CASE(caldyn_thermo) 107 CASE(thermo_theta) 108 theta = temp(ij,l)*(p_ij/preff)**(-kappa) 109 theta_rhodz(ij,l) = mass * theta 110 CASE(thermo_entropy) 111 nu = log(p_ij/preff) 112 chi = log(temp(ij,l)/Treff) 113 entropy = cppd*chi-Rd*nu 114 theta_rhodz(ij,l) = mass * entropy 115 ! CASE(thermo_moist) 116 ! q_ij=q(ij,l,1) 117 ! r_ij=1.-q_ij 118 ! mass=mass*(1-q_ij) ! dry mass 119 ! nu = log(p_ij/preff) 120 ! chi = log(temp(ij,l)/Treff) 121 ! entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu) 122 ! theta_rhodz(ij,l) = mass * entropy 123 CASE DEFAULT 124 STOP 125 END SELECT 126 ENDDO 127 ENDDO 128 ENDDO 129 END SUBROUTINE compute_temperature2entropy 130 131 SUBROUTINE compute_etat0_collocated_hex(ps,mass,phis,temp_i,u, geopot,W, q) 132 USE wind_mod 133 USE disvert_mod 81 SUBROUTINE compute_etat0_collocated_hex(phis,ps,mass,theta_rhodz,u, geopot,W, q) 134 82 REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 135 83 REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm) 84 REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 136 85 REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 137 REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm)138 86 REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 139 87 REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1) … … 141 89 REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) 142 90 143 REAL(rstd) :: ulon_i(iim*jjm,llm)144 REAL(rstd) :: ulat_i(iim*jjm,llm)145 146 91 REAL(rstd) :: ps_e(3*iim*jjm) 92 REAL(rstd) :: phis_e(3*iim*jjm) 147 93 REAL(rstd) :: mass_e(3*iim*jjm,llm) 148 REAL(rstd) :: phis_e(3*iim*jjm)149 REAL(rstd) :: temp_e(3*iim*jjm,llm)150 94 REAL(rstd) :: geopot_e(3*iim*jjm,llm+1) 151 REAL(rstd) :: ulon_e(3*iim*jjm,llm)152 REAL(rstd) :: ulat_e(3*iim*jjm,llm)153 95 REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot) 154 96 155 INTEGER :: l,i,j,ij156 REAL :: p_ik, v_ik, mass_ik157 158 ! For NH geopotential and vertical momentum must be initialized.159 ! Unless autoinit_NH is set to .FALSE. , they will be initialized160 ! to hydrostatic geopotential and zero161 autoinit_mass = .TRUE.162 autoinit_NH = .NOT. hydrostatic163 97 w(:,:) = 0 164 165 CALL compute_etat0_collocated(iim*jjm , lon_i, lat_i, phis, ps, mass, temp_i, ulon_i, ulat_i, geopot, q) 166 CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) 167 168 IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps 169 IF(autoinit_NH) THEN 170 geopot(:,1) = phis(:) ! surface geopotential 171 DO l = 1, llm 172 DO ij=1,iim*jjm 173 ! hybrid pressure coordinate 174 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 175 mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 176 ! v=R.T/p, R=kappa*cpp 177 v_ik = kappa*cpp*temp_i(ij,l)/p_ik 178 geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g 179 END DO 180 END DO 181 END IF 182 183 CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u) 98 CALL compute_etat0_collocated(iim*jjm , lon_i, lat_i, phis, ps, mass, theta_rhodz, geopot, q) 99 CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u) 184 100 185 101 END SUBROUTINE compute_etat0_collocated_hex 186 102 187 SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, temp, ulon, ulat, geopot, q) 103 SUBROUTINE compute_etat0_collocated_unst(phis,ps,mass,theta_rhodz,u, geopot,W, q) 104 REAL(rstd),INTENT(INOUT) :: ps(primal_num) 105 REAL(rstd),INTENT(INOUT) :: mass(llm,primal_num) 106 REAL(rstd),INTENT(OUT) :: theta_rhodz(llm,primal_num) 107 REAL(rstd),INTENT(OUT) :: phis(primal_num) 108 REAL(rstd),INTENT(OUT) :: u(llm, edge_num) 109 REAL(rstd),INTENT(OUT) :: W(llm+1, primal_num) 110 REAL(rstd),INTENT(OUT) :: geopot(llm+1, primal_num) 111 REAL(rstd),INTENT(OUT) :: q(llm, primal_num, nqtot) 112 ! The 2D/3D arrays below have the shape expected by compute_etat0_collocated 113 REAL(rstd) :: ps_e(edge_num) 114 REAL(rstd) :: phis_e(edge_num) 115 REAL(rstd) :: u_e(edge_num, llm) 116 REAL(rstd) :: mass_i(primal_num, llm), theta_rhodz_i(primal_num, llm), mass_e(edge_num, llm) 117 REAL(rstd) :: geopot_i(edge_num, llm+1), geopot_e(edge_num, llm+1) 118 REAL(rstd) :: q_i(primal_num, llm, nqtot), q_e(edge_num, llm, nqtot) 119 INTEGER :: iq 120 121 w(:,:) = 0 122 CALL compute_etat0_collocated(primal_num , lon_i, lat_i, phis, ps, mass_i, theta_rhodz_i, geopot_i, q_i) 123 CALL compute_etat0_collocated(edge_num, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u_e) 124 mass = TRANSPOSE(mass_i) 125 theta_rhodz = TRANSPOSE(theta_rhodz_i) 126 geopot = TRANSPOSE(geopot_i) 127 u = TRANSPOSE(u_e) 128 DO iq=1,nqtot 129 q(:,:,iq) = TRANSPOSE(q_i(:,:,iq)) 130 END DO 131 END SUBROUTINE compute_etat0_collocated_unst 132 133 SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, theta_rhodz, geopot, q, ep, uperp) 188 134 USE etat0_isothermal_mod, ONLY : compute_isothermal => compute_etat0 189 135 USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0 … … 199 145 USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0 200 146 USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0 147 USE disvert_mod, ONLY : ptop, ap, bp, mass_ak, mass_bk, mass_dak, mass_dbk 201 148 INTEGER :: ngrid 202 149 REAL(rstd),INTENT(IN) :: lon(ngrid), lat(ngrid) 203 150 REAL(rstd),INTENT(INOUT) :: ps(ngrid) 204 151 REAL(rstd),INTENT(INOUT) :: mass(ngrid,llm) 152 REAL(rstd),INTENT(OUT) :: theta_rhodz(ngrid,llm) 205 153 REAL(rstd),INTENT(OUT) :: phis(ngrid) 206 REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)207 REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)208 REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)209 154 REAL(rstd),INTENT(OUT) :: geopot(ngrid,llm+1) 210 155 REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) 156 REAL(rstd),INTENT(OUT), OPTIONAL :: ep(ngrid,3), uperp(ngrid,llm) 157 158 REAL(rstd) :: ulon(ngrid,llm) 159 REAL(rstd) :: ulat(ngrid,llm) 160 REAL(rstd) :: temp(ngrid,llm) 161 162 LOGICAL :: autoinit_mass, autoinit_NH 163 INTEGER :: ij,l 164 REAL(rstd) :: clat, slat, clon, slon, u3d(3), p_ik, mass_ik, v_ik 165 REAL(rstd) :: q_ik, r_ik, chi, entropy, theta, log_p_preff 166 167 ! For NH, geopotential and vertical momentum must be initialized. 168 ! Unless autoinit_NH is set to .FALSE. , they will be initialized 169 ! to hydrostatic geopotential and zero 170 autoinit_mass = .TRUE. 171 autoinit_NH = .NOT. hydrostatic 211 172 212 173 SELECT CASE (TRIM(etat0_type)) … … 232 193 ! autoinit_NH = .FALSE. ! compute_bubble initializes geopot 233 194 CASE('williamson91.6') 234 CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), t emp(:,1), ulon(:,1), ulat(:,1))195 CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), theta_rhodz(:,1), ulon(:,1), ulat(:,1)) 235 196 autoinit_mass = .FALSE. ! do not overwrite mass 236 197 CASE('dcmip2016_baroclinic_wave') … … 242 203 END SELECT 243 204 205 IF(PRESENT(uperp)) THEN 206 DO l = 1, llm 207 DO ij=1,ngrid 208 clat = COS(lat(ij)) 209 slat = SIN(lat(ij)) 210 clon = COS(lon(ij)) 211 slon = SIN(lon(ij)) 212 u3d(1) = -clon*slat*ulat(ij,l) - slon*ulon(ij,l) 213 u3d(2) = -slon*slat*ulat(ij,l) + clon*ulon(ij,l) 214 u3d(3) = clat*ulat(ij,l) 215 uperp(ij,l) = SUM(u3d*ep(ij,:)) 216 END DO 217 END DO 218 END IF 219 220 IF(autoinit_mass) THEN 221 DO l = 1, llm 222 DO ij=1,ngrid 223 mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 224 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 225 mass(ij,l) = mass_ik 226 SELECT CASE(caldyn_thermo) 227 CASE(thermo_theta) 228 theta = temp(ij,l)*(p_ik/preff)**(-kappa) 229 theta_rhodz(ij,l) = mass_ik * theta 230 CASE(thermo_entropy) 231 log_p_preff = log(p_ik/preff) 232 chi = log(temp(ij,l)/Treff) 233 entropy = cpp*chi-Rd*log_p_preff 234 theta_rhodz(ij,l) = mass_ik * entropy 235 CASE(thermo_moist) 236 q_ik=q(ij,l,1) 237 r_ik=1.-q_ik 238 mass_ik = mass_ik*(1.-q_ik) ! dry mass 239 log_p_preff = log(p_ik/preff) 240 chi = log(temp(ij,l)/Treff) 241 entropy = r_ik*(cpp*chi-Rd*nu) + q_ik*(cppv*chi-Rv*log_p_preff) 242 theta_rhodz(ij,l) = mass_ik * entropy 243 CASE DEFAULT 244 STOP 245 END SELECT 246 END DO 247 END DO 248 END IF 249 250 IF(autoinit_NH) THEN 251 geopot(:,1) = phis(:) ! surface geopotential 252 DO l = 1, llm 253 DO ij=1,ngrid 254 ! hybrid pressure coordinate 255 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 256 mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 257 ! v=R.T/p, R=kappa*cpp 258 v_ik = kappa*cpp*temp(ij,l)/p_ik 259 geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g 260 END DO 261 END DO 262 END IF 263 244 264 END SUBROUTINE compute_etat0_collocated 245 265
Note: See TracChangeset
for help on using the changeset viewer.
