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