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