| 1 | !IDEAL:MODEL_LAYER:INITIALIZATION |
|---|
| 2 | ! |
|---|
| 3 | |
|---|
| 4 | ! This MODULE holds the routines which are used to perform various initializations |
|---|
| 5 | ! for the individual domains. |
|---|
| 6 | |
|---|
| 7 | ! This MODULE CONTAINS the following routines: |
|---|
| 8 | |
|---|
| 9 | ! initialize_field_test - 1. Set different fields to different constant |
|---|
| 10 | ! values. This is only a test. If the correct |
|---|
| 11 | ! domain is not found (based upon the "id") |
|---|
| 12 | ! then a fatal error is issued. |
|---|
| 13 | |
|---|
| 14 | !----------------------------------------------------------------------- |
|---|
| 15 | |
|---|
| 16 | MODULE module_initialize_ideal |
|---|
| 17 | |
|---|
| 18 | USE module_domain |
|---|
| 19 | USE module_io_domain |
|---|
| 20 | USE module_state_description |
|---|
| 21 | USE module_model_constants |
|---|
| 22 | USE module_bc |
|---|
| 23 | USE module_timing |
|---|
| 24 | USE module_configure |
|---|
| 25 | USE module_init_utilities |
|---|
| 26 | #ifdef DM_PARALLEL |
|---|
| 27 | USE module_dm |
|---|
| 28 | #endif |
|---|
| 29 | |
|---|
| 30 | |
|---|
| 31 | CONTAINS |
|---|
| 32 | |
|---|
| 33 | |
|---|
| 34 | !------------------------------------------------------------------- |
|---|
| 35 | ! this is a wrapper for the solver-specific init_domain routines. |
|---|
| 36 | ! Also dereferences the grid variables and passes them down as arguments. |
|---|
| 37 | ! This is crucial, since the lower level routines may do message passing |
|---|
| 38 | ! and this will get fouled up on machines that insist on passing down |
|---|
| 39 | ! copies of assumed-shape arrays (by passing down as arguments, the |
|---|
| 40 | ! data are treated as assumed-size -- ie. f77 -- arrays and the copying |
|---|
| 41 | ! business is avoided). Fie on the F90 designers. Fie and a pox. |
|---|
| 42 | ! NOTE: Modified to remove all but arrays of rank 4 or more from the |
|---|
| 43 | ! argument list. Arrays with rank>3 are still problematic due to the |
|---|
| 44 | ! above-noted fie- and pox-ities. TBH 20061129. |
|---|
| 45 | |
|---|
| 46 | SUBROUTINE init_domain ( grid ) |
|---|
| 47 | |
|---|
| 48 | IMPLICIT NONE |
|---|
| 49 | |
|---|
| 50 | ! Input data. |
|---|
| 51 | TYPE (domain), POINTER :: grid |
|---|
| 52 | ! Local data. |
|---|
| 53 | INTEGER :: idum1, idum2 |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 ) |
|---|
| 57 | |
|---|
| 58 | CALL init_domain_rk( grid & |
|---|
| 59 | ! |
|---|
| 60 | #include <actual_new_args.inc> |
|---|
| 61 | ! |
|---|
| 62 | ) |
|---|
| 63 | |
|---|
| 64 | END SUBROUTINE init_domain |
|---|
| 65 | |
|---|
| 66 | !------------------------------------------------------------------- |
|---|
| 67 | |
|---|
| 68 | SUBROUTINE init_domain_rk ( grid & |
|---|
| 69 | ! |
|---|
| 70 | # include <dummy_new_args.inc> |
|---|
| 71 | ! |
|---|
| 72 | ) |
|---|
| 73 | IMPLICIT NONE |
|---|
| 74 | |
|---|
| 75 | ! Input data. |
|---|
| 76 | TYPE (domain), POINTER :: grid |
|---|
| 77 | |
|---|
| 78 | # include <dummy_new_decl.inc> |
|---|
| 79 | |
|---|
| 80 | TYPE (grid_config_rec_type) :: config_flags |
|---|
| 81 | |
|---|
| 82 | ! Local data |
|---|
| 83 | INTEGER :: & |
|---|
| 84 | ids, ide, jds, jde, kds, kde, & |
|---|
| 85 | ims, ime, jms, jme, kms, kme, & |
|---|
| 86 | its, ite, jts, jte, kts, kte, & |
|---|
| 87 | i, j, k |
|---|
| 88 | |
|---|
| 89 | ! Local data |
|---|
| 90 | |
|---|
| 91 | INTEGER, PARAMETER :: nl_max = 1000 |
|---|
| 92 | REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in |
|---|
| 93 | INTEGER :: nl_in |
|---|
| 94 | |
|---|
| 95 | |
|---|
| 96 | INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc |
|---|
| 97 | REAL :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u |
|---|
| 98 | REAL :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2 |
|---|
| 99 | ! REAL, EXTERNAL :: interp_0 |
|---|
| 100 | REAL :: hm, xa |
|---|
| 101 | REAL :: pi |
|---|
| 102 | |
|---|
| 103 | ! stuff from original initialization that has been dropped from the Registry |
|---|
| 104 | REAL :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt |
|---|
| 105 | REAL :: qvf1, qvf2, pd_surf |
|---|
| 106 | INTEGER :: it |
|---|
| 107 | real :: thtmp, ptmp, temp(3) |
|---|
| 108 | |
|---|
| 109 | LOGICAL :: moisture_init |
|---|
| 110 | LOGICAL :: stretch_grid, dry_sounding |
|---|
| 111 | |
|---|
| 112 | REAL :: xa1, xal1,pii,hm1 ! data for intercomparison setup from dale |
|---|
| 113 | |
|---|
| 114 | SELECT CASE ( model_data_order ) |
|---|
| 115 | CASE ( DATA_ORDER_ZXY ) |
|---|
| 116 | kds = grid%sd31 ; kde = grid%ed31 ; |
|---|
| 117 | ids = grid%sd32 ; ide = grid%ed32 ; |
|---|
| 118 | jds = grid%sd33 ; jde = grid%ed33 ; |
|---|
| 119 | |
|---|
| 120 | kms = grid%sm31 ; kme = grid%em31 ; |
|---|
| 121 | ims = grid%sm32 ; ime = grid%em32 ; |
|---|
| 122 | jms = grid%sm33 ; jme = grid%em33 ; |
|---|
| 123 | |
|---|
| 124 | kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch |
|---|
| 125 | its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch |
|---|
| 126 | jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch |
|---|
| 127 | CASE ( DATA_ORDER_XYZ ) |
|---|
| 128 | ids = grid%sd31 ; ide = grid%ed31 ; |
|---|
| 129 | jds = grid%sd32 ; jde = grid%ed32 ; |
|---|
| 130 | kds = grid%sd33 ; kde = grid%ed33 ; |
|---|
| 131 | |
|---|
| 132 | ims = grid%sm31 ; ime = grid%em31 ; |
|---|
| 133 | jms = grid%sm32 ; jme = grid%em32 ; |
|---|
| 134 | kms = grid%sm33 ; kme = grid%em33 ; |
|---|
| 135 | |
|---|
| 136 | its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch |
|---|
| 137 | jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch |
|---|
| 138 | kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch |
|---|
| 139 | CASE ( DATA_ORDER_XZY ) |
|---|
| 140 | ids = grid%sd31 ; ide = grid%ed31 ; |
|---|
| 141 | kds = grid%sd32 ; kde = grid%ed32 ; |
|---|
| 142 | jds = grid%sd33 ; jde = grid%ed33 ; |
|---|
| 143 | |
|---|
| 144 | ims = grid%sm31 ; ime = grid%em31 ; |
|---|
| 145 | kms = grid%sm32 ; kme = grid%em32 ; |
|---|
| 146 | jms = grid%sm33 ; jme = grid%em33 ; |
|---|
| 147 | |
|---|
| 148 | its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch |
|---|
| 149 | kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch |
|---|
| 150 | jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch |
|---|
| 151 | |
|---|
| 152 | END SELECT |
|---|
| 153 | |
|---|
| 154 | |
|---|
| 155 | hm = 100. |
|---|
| 156 | xa = 5.0 |
|---|
| 157 | |
|---|
| 158 | icm = ide/2 |
|---|
| 159 | |
|---|
| 160 | |
|---|
| 161 | xa1 = 5000./500. |
|---|
| 162 | xal1 = 4000./500. |
|---|
| 163 | pii = 2.*asin(1.0) |
|---|
| 164 | hm1 = 250. |
|---|
| 165 | ! hm1 = 1000. |
|---|
| 166 | |
|---|
| 167 | |
|---|
| 168 | stretch_grid = .true. |
|---|
| 169 | ! z_scale = .50 |
|---|
| 170 | z_scale = .40 |
|---|
| 171 | pi = 2.*asin(1.0) |
|---|
| 172 | write(6,*) ' pi is ',pi |
|---|
| 173 | nxc = (ide-ids)/4 |
|---|
| 174 | nyc = (jde-jds)/2 |
|---|
| 175 | |
|---|
| 176 | CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) |
|---|
| 177 | |
|---|
| 178 | ! here we check to see if the boundary conditions are set properly |
|---|
| 179 | |
|---|
| 180 | CALL boundary_condition_check( config_flags, bdyzone, error, grid%id ) |
|---|
| 181 | |
|---|
| 182 | moisture_init = .true. |
|---|
| 183 | |
|---|
| 184 | grid%itimestep=0 |
|---|
| 185 | |
|---|
| 186 | #ifdef DM_PARALLEL |
|---|
| 187 | CALL wrf_dm_bcast_bytes( icm , IWORDSIZE ) |
|---|
| 188 | CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE ) |
|---|
| 189 | #endif |
|---|
| 190 | |
|---|
| 191 | CALL nl_set_mminlu(1,' ') |
|---|
| 192 | CALL nl_set_iswater(1,0) |
|---|
| 193 | CALL nl_set_cen_lat(1,40.) |
|---|
| 194 | CALL nl_set_cen_lon(1,-105.) |
|---|
| 195 | CALL nl_set_truelat1(1,0.) |
|---|
| 196 | CALL nl_set_truelat2(1,0.) |
|---|
| 197 | CALL nl_set_moad_cen_lat (1,0.) |
|---|
| 198 | CALL nl_set_stand_lon (1,0.) |
|---|
| 199 | CALL nl_set_map_proj(1,0) |
|---|
| 200 | |
|---|
| 201 | |
|---|
| 202 | ! here we initialize data we currently is not initialized |
|---|
| 203 | ! in the input data |
|---|
| 204 | |
|---|
| 205 | DO j = jts, jte |
|---|
| 206 | DO i = its, ite |
|---|
| 207 | grid%msftx(i,j) = 1. |
|---|
| 208 | grid%msfty(i,j) = 1. |
|---|
| 209 | grid%msfux(i,j) = 1. |
|---|
| 210 | grid%msfuy(i,j) = 1. |
|---|
| 211 | grid%msfvx(i,j) = 1. |
|---|
| 212 | grid%msfvx_inv(i,j)= 1. |
|---|
| 213 | grid%msfvy(i,j) = 1. |
|---|
| 214 | grid%sina(i,j) = 0. |
|---|
| 215 | grid%cosa(i,j) = 1. |
|---|
| 216 | grid%e(i,j) = 0. |
|---|
| 217 | grid%f(i,j) = 0. |
|---|
| 218 | |
|---|
| 219 | END DO |
|---|
| 220 | END DO |
|---|
| 221 | |
|---|
| 222 | DO j = jts, jte |
|---|
| 223 | DO k = kts, kte |
|---|
| 224 | DO i = its, ite |
|---|
| 225 | grid%ww(i,k,j) = 0. |
|---|
| 226 | END DO |
|---|
| 227 | END DO |
|---|
| 228 | END DO |
|---|
| 229 | |
|---|
| 230 | grid%step_number = 0 |
|---|
| 231 | |
|---|
| 232 | ! set up the grid |
|---|
| 233 | |
|---|
| 234 | IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz) |
|---|
| 235 | DO k=1, kde |
|---|
| 236 | grid%znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ & |
|---|
| 237 | (1.-exp(-1./z_scale)) |
|---|
| 238 | ENDDO |
|---|
| 239 | ELSE |
|---|
| 240 | DO k=1, kde |
|---|
| 241 | grid%znw(k) = 1. - float(k-1)/float(kde-1) |
|---|
| 242 | ENDDO |
|---|
| 243 | ENDIF |
|---|
| 244 | |
|---|
| 245 | DO k=1, kde-1 |
|---|
| 246 | grid%dnw(k) = grid%znw(k+1) - grid%znw(k) |
|---|
| 247 | grid%rdnw(k) = 1./grid%dnw(k) |
|---|
| 248 | grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k)) |
|---|
| 249 | ENDDO |
|---|
| 250 | DO k=2, kde-1 |
|---|
| 251 | grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1)) |
|---|
| 252 | grid%rdn(k) = 1./grid%dn(k) |
|---|
| 253 | grid%fnp(k) = .5* grid%dnw(k )/grid%dn(k) |
|---|
| 254 | grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k) |
|---|
| 255 | ENDDO |
|---|
| 256 | |
|---|
| 257 | cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2) |
|---|
| 258 | cof2 = grid%dn(2) /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3) |
|---|
| 259 | grid%cf1 = grid%fnp(2) + cof1 |
|---|
| 260 | grid%cf2 = grid%fnm(2) - cof1 - cof2 |
|---|
| 261 | grid%cf3 = cof2 |
|---|
| 262 | |
|---|
| 263 | grid%cfn = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1) |
|---|
| 264 | grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1) |
|---|
| 265 | grid%rdx = 1./config_flags%dx |
|---|
| 266 | grid%rdy = 1./config_flags%dy |
|---|
| 267 | |
|---|
| 268 | ! get the sounding from the ascii sounding file, first get dry sounding and |
|---|
| 269 | ! calculate base state |
|---|
| 270 | |
|---|
| 271 | write(6,*) ' getting dry sounding for base state ' |
|---|
| 272 | dry_sounding = .true. |
|---|
| 273 | CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, & |
|---|
| 274 | nl_max, nl_in, .true.) |
|---|
| 275 | |
|---|
| 276 | write(6,*) ' returned from reading sounding, nl_in is ',nl_in |
|---|
| 277 | |
|---|
| 278 | |
|---|
| 279 | ! find ptop for the desired ztop (ztop is input from the namelist), |
|---|
| 280 | ! and find surface pressure |
|---|
| 281 | |
|---|
| 282 | grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in ) |
|---|
| 283 | |
|---|
| 284 | DO j=jts,jte |
|---|
| 285 | DO i=its,ite ! flat surface |
|---|
| 286 | !! grid%ht(i,j) = 0. |
|---|
| 287 | grid%ht(i,j) = hm/(1.+(float(i-icm)/xa)**2) |
|---|
| 288 | ! grid%ht(i,j) = hm1*exp(-(( float(i-icm)/xa1)**2)) & |
|---|
| 289 | ! *( (cos(pii*float(i-icm)/xal1))**2 ) |
|---|
| 290 | grid%phb(i,1,j) = g*grid%ht(i,j) |
|---|
| 291 | grid%php(i,1,j) = 0. |
|---|
| 292 | grid%ph0(i,1,j) = grid%phb(i,1,j) |
|---|
| 293 | ENDDO |
|---|
| 294 | ENDDO |
|---|
| 295 | |
|---|
| 296 | DO J = jts, jte |
|---|
| 297 | DO I = its, ite |
|---|
| 298 | |
|---|
| 299 | p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in ) |
|---|
| 300 | grid%mub(i,j) = p_surf-grid%p_top |
|---|
| 301 | |
|---|
| 302 | ! this is dry hydrostatic sounding (base state), so given grid%p (coordinate), |
|---|
| 303 | ! interp theta (from interp) and compute 1/rho from eqn. of state |
|---|
| 304 | |
|---|
| 305 | DO K = 1, kte-1 |
|---|
| 306 | p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top |
|---|
| 307 | grid%pb(i,k,j) = p_level |
|---|
| 308 | grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0 |
|---|
| 309 | grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm |
|---|
| 310 | ENDDO |
|---|
| 311 | |
|---|
| 312 | ! calc hydrostatic balance (alternatively we could interp the geopotential from the |
|---|
| 313 | ! sounding, but this assures that the base state is in exact hydrostatic balance with |
|---|
| 314 | ! respect to the model eqns. |
|---|
| 315 | |
|---|
| 316 | DO k = 2,kte |
|---|
| 317 | grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j) |
|---|
| 318 | ENDDO |
|---|
| 319 | |
|---|
| 320 | ENDDO |
|---|
| 321 | ENDDO |
|---|
| 322 | |
|---|
| 323 | write(6,*) ' ptop is ',grid%p_top |
|---|
| 324 | write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top |
|---|
| 325 | |
|---|
| 326 | ! calculate full state for each column - this includes moisture. |
|---|
| 327 | |
|---|
| 328 | write(6,*) ' getting moist sounding for full state ' |
|---|
| 329 | dry_sounding = .false. |
|---|
| 330 | CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, & |
|---|
| 331 | nl_max, nl_in, .false. ) |
|---|
| 332 | |
|---|
| 333 | DO J = jts, min(jde-1,jte) |
|---|
| 334 | DO I = its, min(ide-1,ite) |
|---|
| 335 | |
|---|
| 336 | ! At this point grid%p_top is already set. find the DRY mass in the column |
|---|
| 337 | ! by interpolating the DRY pressure. |
|---|
| 338 | |
|---|
| 339 | pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in ) |
|---|
| 340 | |
|---|
| 341 | ! compute the perturbation mass and the full mass |
|---|
| 342 | |
|---|
| 343 | grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j) |
|---|
| 344 | grid%mu_2(i,j) = grid%mu_1(i,j) |
|---|
| 345 | grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j) |
|---|
| 346 | |
|---|
| 347 | ! given the dry pressure and coordinate system, interp the potential |
|---|
| 348 | ! temperature and qv |
|---|
| 349 | |
|---|
| 350 | do k=1,kde-1 |
|---|
| 351 | |
|---|
| 352 | p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top |
|---|
| 353 | |
|---|
| 354 | moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in ) |
|---|
| 355 | grid%t_1(i,k,j) = interp_0( theta, pd_in, p_level, nl_in ) - t0 |
|---|
| 356 | grid%t_2(i,k,j) = grid%t_1(i,k,j) |
|---|
| 357 | |
|---|
| 358 | |
|---|
| 359 | enddo |
|---|
| 360 | |
|---|
| 361 | ! integrate the hydrostatic equation (from the RHS of the bigstep |
|---|
| 362 | ! vertical momentum equation) down from the top to get grid%p. |
|---|
| 363 | ! first from the top of the model to the top pressure |
|---|
| 364 | |
|---|
| 365 | k = kte-1 ! top level |
|---|
| 366 | |
|---|
| 367 | qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV)) |
|---|
| 368 | qvf2 = 1./(1.+qvf1) |
|---|
| 369 | qvf1 = qvf1*qvf2 |
|---|
| 370 | |
|---|
| 371 | ! grid%p(i,k,j) = - 0.5*grid%mu_1(i,j)/grid%rdnw(k) |
|---|
| 372 | grid%p(i,k,j) = - 0.5*(grid%mu_1(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2 |
|---|
| 373 | qvf = 1. + rvovrd*moist(i,k,j,P_QV) |
|---|
| 374 | grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* & |
|---|
| 375 | (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm) |
|---|
| 376 | grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j) |
|---|
| 377 | |
|---|
| 378 | ! down the column |
|---|
| 379 | |
|---|
| 380 | do k=kte-2,1,-1 |
|---|
| 381 | qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV)) |
|---|
| 382 | qvf2 = 1./(1.+qvf1) |
|---|
| 383 | qvf1 = qvf1*qvf2 |
|---|
| 384 | grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_1(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1) |
|---|
| 385 | qvf = 1. + rvovrd*moist(i,k,j,P_QV) |
|---|
| 386 | grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* & |
|---|
| 387 | (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm) |
|---|
| 388 | grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j) |
|---|
| 389 | enddo |
|---|
| 390 | |
|---|
| 391 | ! this is the hydrostatic equation used in the model after the |
|---|
| 392 | ! small timesteps. In the model, grid%al (inverse density) |
|---|
| 393 | ! is computed from the geopotential. |
|---|
| 394 | |
|---|
| 395 | |
|---|
| 396 | grid%ph_1(i,1,j) = 0. |
|---|
| 397 | DO k = 2,kte |
|---|
| 398 | grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*( & |
|---|
| 399 | (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ & |
|---|
| 400 | grid%mu_1(i,j)*grid%alb(i,k-1,j) ) |
|---|
| 401 | |
|---|
| 402 | grid%ph_2(i,k,j) = grid%ph_1(i,k,j) |
|---|
| 403 | grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) |
|---|
| 404 | ENDDO |
|---|
| 405 | |
|---|
| 406 | if((i==2) .and. (j==2)) then |
|---|
| 407 | write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),& |
|---|
| 408 | grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), & |
|---|
| 409 | grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1) |
|---|
| 410 | endif |
|---|
| 411 | |
|---|
| 412 | ENDDO |
|---|
| 413 | ENDDO |
|---|
| 414 | |
|---|
| 415 | write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1) |
|---|
| 416 | write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv ' |
|---|
| 417 | do k=1,kde-1 |
|---|
| 418 | write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), & |
|---|
| 419 | grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), & |
|---|
| 420 | grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV) |
|---|
| 421 | enddo |
|---|
| 422 | |
|---|
| 423 | write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv ' |
|---|
| 424 | do k=1,kde-1 |
|---|
| 425 | write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), & |
|---|
| 426 | grid%p(1,k,1), grid%al(1,k,1), & |
|---|
| 427 | grid%t_1(1,k,1), moist(1,k,1,P_QV) |
|---|
| 428 | enddo |
|---|
| 429 | |
|---|
| 430 | ! interp v |
|---|
| 431 | |
|---|
| 432 | DO J = jts, jte |
|---|
| 433 | DO I = its, min(ide-1,ite) |
|---|
| 434 | |
|---|
| 435 | IF (j == jds) THEN |
|---|
| 436 | z_at_v = grid%phb(i,1,j)/g |
|---|
| 437 | ELSE IF (j == jde) THEN |
|---|
| 438 | z_at_v = grid%phb(i,1,j-1)/g |
|---|
| 439 | ELSE |
|---|
| 440 | z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g |
|---|
| 441 | END IF |
|---|
| 442 | |
|---|
| 443 | p_surf = interp_0( p_in, zk, z_at_v, nl_in ) |
|---|
| 444 | |
|---|
| 445 | DO K = 1, kte |
|---|
| 446 | p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top |
|---|
| 447 | grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in ) |
|---|
| 448 | grid%v_2(i,k,j) = grid%v_1(i,k,j) |
|---|
| 449 | ENDDO |
|---|
| 450 | |
|---|
| 451 | ENDDO |
|---|
| 452 | ENDDO |
|---|
| 453 | |
|---|
| 454 | ! interp u |
|---|
| 455 | |
|---|
| 456 | DO J = jts, min(jde-1,jte) |
|---|
| 457 | DO I = its, ite |
|---|
| 458 | |
|---|
| 459 | IF (i == ids) THEN |
|---|
| 460 | z_at_u = grid%phb(i,1,j)/g |
|---|
| 461 | ELSE IF (i == ide) THEN |
|---|
| 462 | z_at_u = grid%phb(i-1,1,j)/g |
|---|
| 463 | ELSE |
|---|
| 464 | z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g |
|---|
| 465 | END IF |
|---|
| 466 | |
|---|
| 467 | p_surf = interp_0( p_in, zk, z_at_u, nl_in ) |
|---|
| 468 | |
|---|
| 469 | DO K = 1, kte |
|---|
| 470 | p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top |
|---|
| 471 | grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in ) |
|---|
| 472 | grid%u_2(i,k,j) = grid%u_1(i,k,j) |
|---|
| 473 | ENDDO |
|---|
| 474 | |
|---|
| 475 | ENDDO |
|---|
| 476 | ENDDO |
|---|
| 477 | |
|---|
| 478 | ! set w |
|---|
| 479 | |
|---|
| 480 | DO J = jts, min(jde-1,jte) |
|---|
| 481 | DO K = kts, kte |
|---|
| 482 | DO I = its, min(ide-1,ite) |
|---|
| 483 | grid%w_1(i,k,j) = 0. |
|---|
| 484 | grid%w_2(i,k,j) = 0. |
|---|
| 485 | ENDDO |
|---|
| 486 | ENDDO |
|---|
| 487 | ENDDO |
|---|
| 488 | |
|---|
| 489 | ! set a few more things |
|---|
| 490 | |
|---|
| 491 | DO J = jts, min(jde-1,jte) |
|---|
| 492 | DO K = kts, kte-1 |
|---|
| 493 | DO I = its, min(ide-1,ite) |
|---|
| 494 | grid%h_diabatic(i,k,j) = 0. |
|---|
| 495 | ENDDO |
|---|
| 496 | ENDDO |
|---|
| 497 | ENDDO |
|---|
| 498 | |
|---|
| 499 | DO k=1,kte-1 |
|---|
| 500 | grid%t_base(k) = grid%t_1(1,k,1) |
|---|
| 501 | grid%qv_base(k) = moist(1,k,1,P_QV) |
|---|
| 502 | grid%u_base(k) = grid%u_1(1,k,1) |
|---|
| 503 | grid%v_base(k) = grid%v_1(1,k,1) |
|---|
| 504 | grid%z_base(k) = 0.5*(grid%phb(1,k,1)+grid%phb(1,k+1,1)+grid%ph_1(1,k,1)+grid%ph_1(1,k+1,1))/g |
|---|
| 505 | ENDDO |
|---|
| 506 | |
|---|
| 507 | DO J = jts, min(jde-1,jte) |
|---|
| 508 | DO I = its, min(ide-1,ite) |
|---|
| 509 | thtmp = grid%t_2(i,1,j)+t0 |
|---|
| 510 | ptmp = grid%p(i,1,j)+grid%pb(i,1,j) |
|---|
| 511 | temp(1) = thtmp * (ptmp/p1000mb)**rcp |
|---|
| 512 | thtmp = grid%t_2(i,2,j)+t0 |
|---|
| 513 | ptmp = grid%p(i,2,j)+grid%pb(i,2,j) |
|---|
| 514 | temp(2) = thtmp * (ptmp/p1000mb)**rcp |
|---|
| 515 | thtmp = grid%t_2(i,3,j)+t0 |
|---|
| 516 | ptmp = grid%p(i,3,j)+grid%pb(i,3,j) |
|---|
| 517 | temp(3) = thtmp * (ptmp/p1000mb)**rcp |
|---|
| 518 | |
|---|
| 519 | grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3) |
|---|
| 520 | grid%tmn(I,J)=grid%tsk(I,J)-0.5 |
|---|
| 521 | ENDDO |
|---|
| 522 | ENDDO |
|---|
| 523 | |
|---|
| 524 | RETURN |
|---|
| 525 | |
|---|
| 526 | END SUBROUTINE init_domain_rk |
|---|
| 527 | |
|---|
| 528 | SUBROUTINE init_module_initialize |
|---|
| 529 | END SUBROUTINE init_module_initialize |
|---|
| 530 | |
|---|
| 531 | !--------------------------------------------------------------------- |
|---|
| 532 | |
|---|
| 533 | ! test driver for get_sounding |
|---|
| 534 | ! |
|---|
| 535 | ! implicit none |
|---|
| 536 | ! integer n |
|---|
| 537 | ! parameter(n = 1000) |
|---|
| 538 | ! real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n) |
|---|
| 539 | ! logical dry |
|---|
| 540 | ! integer nl,k |
|---|
| 541 | ! |
|---|
| 542 | ! dry = .false. |
|---|
| 543 | ! dry = .true. |
|---|
| 544 | ! call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl ) |
|---|
| 545 | ! write(6,*) ' input levels ',nl |
|---|
| 546 | ! write(6,*) ' sounding ' |
|---|
| 547 | ! write(6,*) ' k height(m) press (Pa) pd(Pa) theta (K) den(kg/m^3) u(m/s) v(m/s) qv(g/g) ' |
|---|
| 548 | ! do k=1,nl |
|---|
| 549 | ! write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), pd(k), theta(k), rho(k), u(k), v(k), qv(k) |
|---|
| 550 | ! enddo |
|---|
| 551 | ! end |
|---|
| 552 | ! |
|---|
| 553 | !--------------------------------------------------------------------------- |
|---|
| 554 | |
|---|
| 555 | subroutine get_sounding( zk, p, p_dry, theta, rho, & |
|---|
| 556 | u, v, qv, dry, nl_max, nl_in, base_state ) |
|---|
| 557 | implicit none |
|---|
| 558 | |
|---|
| 559 | integer nl_max, nl_in |
|---|
| 560 | real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), & |
|---|
| 561 | u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max) |
|---|
| 562 | logical dry |
|---|
| 563 | logical base_state |
|---|
| 564 | |
|---|
| 565 | integer n, iz |
|---|
| 566 | parameter(n=1000) |
|---|
| 567 | logical debug |
|---|
| 568 | parameter( debug = .false.) |
|---|
| 569 | |
|---|
| 570 | ! input sounding data |
|---|
| 571 | |
|---|
| 572 | real p_surf, th_surf, qv_surf |
|---|
| 573 | real pi_surf, pi(n) |
|---|
| 574 | real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n) |
|---|
| 575 | |
|---|
| 576 | ! diagnostics |
|---|
| 577 | |
|---|
| 578 | real rho_surf, p_input(n), rho_input(n) |
|---|
| 579 | real pm_input(n) ! this are for full moist sounding |
|---|
| 580 | |
|---|
| 581 | ! local data |
|---|
| 582 | |
|---|
| 583 | real p1000mb,cv,cp,r,cvpm,g |
|---|
| 584 | parameter (p1000mb = 1.e+05, r = 287, cp = 1003., cv = cp-r, cvpm = -cv/cp, g=9.81 ) |
|---|
| 585 | integer k, it, nl |
|---|
| 586 | real qvf, qvf1, dz |
|---|
| 587 | |
|---|
| 588 | ! first, read the sounding |
|---|
| 589 | |
|---|
| 590 | call read_sounding( p_surf, th_surf, qv_surf, & |
|---|
| 591 | h_input, th_input, qv_input, u_input, v_input,n, nl, debug ) |
|---|
| 592 | |
|---|
| 593 | ! iz = 1 |
|---|
| 594 | ! do k=2,nl |
|---|
| 595 | ! if(h_input(k) .lt. 12000.) iz = k |
|---|
| 596 | ! enddo |
|---|
| 597 | ! write(6,*) " tropopause ",iz,h_input(iz) |
|---|
| 598 | ! if(dry) then |
|---|
| 599 | ! write(6,*) ' nl is ',nl |
|---|
| 600 | ! do k=1,nl |
|---|
| 601 | ! th_input(k) = th_input(k)+10.+10*float(k)/nl |
|---|
| 602 | ! enddo |
|---|
| 603 | ! write(6,*) ' finished adjusting theta ' |
|---|
| 604 | ! endif |
|---|
| 605 | |
|---|
| 606 | ! do k=1,nl |
|---|
| 607 | ! u_input(k) = 2*u_input(k) |
|---|
| 608 | ! enddo |
|---|
| 609 | ! |
|---|
| 610 | ! end if |
|---|
| 611 | |
|---|
| 612 | if(dry) then |
|---|
| 613 | do k=1,nl |
|---|
| 614 | qv_input(k) = 0. |
|---|
| 615 | enddo |
|---|
| 616 | endif |
|---|
| 617 | |
|---|
| 618 | if(debug) write(6,*) ' number of input levels = ',nl |
|---|
| 619 | |
|---|
| 620 | nl_in = nl |
|---|
| 621 | if(nl_in .gt. nl_max ) then |
|---|
| 622 | write(6,*) ' too many levels for input arrays ',nl_in,nl_max |
|---|
| 623 | call wrf_error_fatal ( ' too many levels for input arrays ' ) |
|---|
| 624 | end if |
|---|
| 625 | |
|---|
| 626 | ! compute diagnostics, |
|---|
| 627 | ! first, convert qv(g/kg) to qv(g/g) |
|---|
| 628 | |
|---|
| 629 | do k=1,nl |
|---|
| 630 | qv_input(k) = 0.001*qv_input(k) |
|---|
| 631 | enddo |
|---|
| 632 | |
|---|
| 633 | p_surf = 100.*p_surf ! convert to pascals |
|---|
| 634 | qvf = 1. + rvovrd*qv_input(1) |
|---|
| 635 | rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm)) |
|---|
| 636 | pi_surf = (p_surf/p1000mb)**(r/cp) |
|---|
| 637 | |
|---|
| 638 | if(debug) then |
|---|
| 639 | write(6,*) ' surface density is ',rho_surf |
|---|
| 640 | write(6,*) ' surface pi is ',pi_surf |
|---|
| 641 | end if |
|---|
| 642 | |
|---|
| 643 | |
|---|
| 644 | ! integrate moist sounding hydrostatically, starting from the |
|---|
| 645 | ! specified surface pressure |
|---|
| 646 | ! -> first, integrate from surface to lowest level |
|---|
| 647 | |
|---|
| 648 | qvf = 1. + rvovrd*qv_input(1) |
|---|
| 649 | qvf1 = 1. + qv_input(1) |
|---|
| 650 | rho_input(1) = rho_surf |
|---|
| 651 | dz = h_input(1) |
|---|
| 652 | do it=1,10 |
|---|
| 653 | pm_input(1) = p_surf & |
|---|
| 654 | - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1 |
|---|
| 655 | rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm)) |
|---|
| 656 | enddo |
|---|
| 657 | |
|---|
| 658 | ! integrate up the column |
|---|
| 659 | |
|---|
| 660 | do k=2,nl |
|---|
| 661 | rho_input(k) = rho_input(k-1) |
|---|
| 662 | dz = h_input(k)-h_input(k-1) |
|---|
| 663 | qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k))) |
|---|
| 664 | qvf = 1. + rvovrd*qv_input(k) ! qv is in g/kg here |
|---|
| 665 | |
|---|
| 666 | do it=1,20 |
|---|
| 667 | pm_input(k) = pm_input(k-1) & |
|---|
| 668 | - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1 |
|---|
| 669 | rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm)) |
|---|
| 670 | enddo |
|---|
| 671 | enddo |
|---|
| 672 | |
|---|
| 673 | ! we have the moist sounding |
|---|
| 674 | |
|---|
| 675 | ! next, compute the dry sounding using p at the highest level from the |
|---|
| 676 | ! moist sounding and integrating down. |
|---|
| 677 | |
|---|
| 678 | p_input(nl) = pm_input(nl) |
|---|
| 679 | |
|---|
| 680 | do k=nl-1,1,-1 |
|---|
| 681 | dz = h_input(k+1)-h_input(k) |
|---|
| 682 | p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g |
|---|
| 683 | enddo |
|---|
| 684 | |
|---|
| 685 | ! write(6,*) ' zeroing u input ' |
|---|
| 686 | |
|---|
| 687 | do k=1,nl |
|---|
| 688 | |
|---|
| 689 | zk(k) = h_input(k) |
|---|
| 690 | p(k) = pm_input(k) |
|---|
| 691 | p_dry(k) = p_input(k) |
|---|
| 692 | theta(k) = th_input(k) |
|---|
| 693 | rho(k) = rho_input(k) |
|---|
| 694 | u(k) = u_input(k) |
|---|
| 695 | ! u(k) = 0. |
|---|
| 696 | v(k) = v_input(k) |
|---|
| 697 | qv(k) = qv_input(k) |
|---|
| 698 | |
|---|
| 699 | enddo |
|---|
| 700 | |
|---|
| 701 | if(debug) then |
|---|
| 702 | write(6,*) ' sounding ' |
|---|
| 703 | write(6,*) ' k height(m) press (Pa) pd(Pa) theta (K) den(kg/m^3) u(m/s) v(m/s) qv(g/g) ' |
|---|
| 704 | do k=1,nl |
|---|
| 705 | write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), p_dry(k), theta(k), rho(k), u(k), v(k), qv(k) |
|---|
| 706 | enddo |
|---|
| 707 | |
|---|
| 708 | end if |
|---|
| 709 | |
|---|
| 710 | end subroutine get_sounding |
|---|
| 711 | |
|---|
| 712 | !------------------------------------------------------- |
|---|
| 713 | |
|---|
| 714 | subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug ) |
|---|
| 715 | implicit none |
|---|
| 716 | integer n,nl |
|---|
| 717 | real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n) |
|---|
| 718 | logical end_of_file |
|---|
| 719 | logical debug |
|---|
| 720 | |
|---|
| 721 | integer k |
|---|
| 722 | |
|---|
| 723 | open(unit=10,file='input_sounding',form='formatted',status='old') |
|---|
| 724 | rewind(10) |
|---|
| 725 | read(10,*) ps, ts, qvs |
|---|
| 726 | if(debug) then |
|---|
| 727 | write(6,*) ' input sounding surface parameters ' |
|---|
| 728 | write(6,*) ' surface pressure (mb) ',ps |
|---|
| 729 | write(6,*) ' surface pot. temp (K) ',ts |
|---|
| 730 | write(6,*) ' surface mixing ratio (g/kg) ',qvs |
|---|
| 731 | end if |
|---|
| 732 | |
|---|
| 733 | end_of_file = .false. |
|---|
| 734 | k = 0 |
|---|
| 735 | |
|---|
| 736 | do while (.not. end_of_file) |
|---|
| 737 | |
|---|
| 738 | read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1) |
|---|
| 739 | k = k+1 |
|---|
| 740 | if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k) |
|---|
| 741 | go to 110 |
|---|
| 742 | 100 end_of_file = .true. |
|---|
| 743 | 110 continue |
|---|
| 744 | enddo |
|---|
| 745 | |
|---|
| 746 | nl = k |
|---|
| 747 | |
|---|
| 748 | close(unit=10,status = 'keep') |
|---|
| 749 | |
|---|
| 750 | end subroutine read_sounding |
|---|
| 751 | |
|---|
| 752 | END MODULE module_initialize_ideal |
|---|
