| 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 |
|---|
| 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 :: dyn_opt |
|---|
| 51 | INTEGER :: idum1, idum2 |
|---|
| 52 | |
|---|
| 53 | CALL nl_get_dyn_opt( 1,dyn_opt ) |
|---|
| 54 | |
|---|
| 55 | CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 ) |
|---|
| 56 | |
|---|
| 57 | IF ( dyn_opt .eq. 1 & |
|---|
| 58 | .or. dyn_opt .eq. 2 & |
|---|
| 59 | .or. dyn_opt .eq. 3 & |
|---|
| 60 | ) THEN |
|---|
| 61 | CALL init_domain_rk( grid & |
|---|
| 62 | ! |
|---|
| 63 | #include <em_actual_new_args.inc> |
|---|
| 64 | ! |
|---|
| 65 | ) |
|---|
| 66 | |
|---|
| 67 | ELSE |
|---|
| 68 | WRITE(0,*)' init_domain: unknown or unimplemented dyn_opt = ',dyn_opt |
|---|
| 69 | CALL wrf_error_fatal ( ' init_domain: unknown or unimplemented dyn_opt ' ) |
|---|
| 70 | ENDIF |
|---|
| 71 | |
|---|
| 72 | END SUBROUTINE init_domain |
|---|
| 73 | |
|---|
| 74 | !------------------------------------------------------------------- |
|---|
| 75 | |
|---|
| 76 | SUBROUTINE init_domain_rk ( grid & |
|---|
| 77 | ! |
|---|
| 78 | # include <em_dummy_new_args.inc> |
|---|
| 79 | ! |
|---|
| 80 | ) |
|---|
| 81 | IMPLICIT NONE |
|---|
| 82 | |
|---|
| 83 | ! Input data. |
|---|
| 84 | TYPE (domain), POINTER :: grid |
|---|
| 85 | |
|---|
| 86 | # include <em_dummy_new_decl.inc> |
|---|
| 87 | |
|---|
| 88 | TYPE (grid_config_rec_type) :: config_flags |
|---|
| 89 | |
|---|
| 90 | ! Local data |
|---|
| 91 | INTEGER :: & |
|---|
| 92 | ids, ide, jds, jde, kds, kde, & |
|---|
| 93 | ims, ime, jms, jme, kms, kme, & |
|---|
| 94 | its, ite, jts, jte, kts, kte, & |
|---|
| 95 | i, j, k |
|---|
| 96 | |
|---|
| 97 | ! Local data |
|---|
| 98 | !****Mars |
|---|
| 99 | REAL :: x_param,y_param,rho_param,dilat |
|---|
| 100 | REAL :: mulu, mulv, addu, addv |
|---|
| 101 | !****Mars |
|---|
| 102 | INTEGER, PARAMETER :: nl_max = 1000 |
|---|
| 103 | REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in |
|---|
| 104 | INTEGER :: nl_in |
|---|
| 105 | |
|---|
| 106 | |
|---|
| 107 | INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc |
|---|
| 108 | REAL :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u |
|---|
| 109 | REAL :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2 |
|---|
| 110 | ! REAL, EXTERNAL :: interp_0 |
|---|
| 111 | REAL :: hm, xa |
|---|
| 112 | REAL :: pi |
|---|
| 113 | |
|---|
| 114 | ! stuff from original initialization that has been dropped from the Registry |
|---|
| 115 | REAL :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt |
|---|
| 116 | REAL :: qvf1, qvf2, pd_surf |
|---|
| 117 | INTEGER :: it |
|---|
| 118 | real :: thtmp, ptmp, temp(3) |
|---|
| 119 | |
|---|
| 120 | LOGICAL :: moisture_init |
|---|
| 121 | LOGICAL :: stretch_grid, dry_sounding |
|---|
| 122 | |
|---|
| 123 | INTEGER :: xs , xe , ys , ye |
|---|
| 124 | REAL :: mtn_ht |
|---|
| 125 | LOGICAL, EXTERNAL :: wrf_dm_on_monitor |
|---|
| 126 | |
|---|
| 127 | !!MARS |
|---|
| 128 | REAL :: lon_input, lat_input, alt_input, tsurf_input |
|---|
| 129 | !!MARS |
|---|
| 130 | |
|---|
| 131 | |
|---|
| 132 | #ifdef DM_PARALLEL |
|---|
| 133 | # include <em_data_calls.inc> |
|---|
| 134 | #endif |
|---|
| 135 | |
|---|
| 136 | |
|---|
| 137 | SELECT CASE ( model_data_order ) |
|---|
| 138 | CASE ( DATA_ORDER_ZXY ) |
|---|
| 139 | kds = grid%sd31 ; kde = grid%ed31 ; |
|---|
| 140 | ids = grid%sd32 ; ide = grid%ed32 ; |
|---|
| 141 | jds = grid%sd33 ; jde = grid%ed33 ; |
|---|
| 142 | |
|---|
| 143 | kms = grid%sm31 ; kme = grid%em31 ; |
|---|
| 144 | ims = grid%sm32 ; ime = grid%em32 ; |
|---|
| 145 | jms = grid%sm33 ; jme = grid%em33 ; |
|---|
| 146 | |
|---|
| 147 | kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch |
|---|
| 148 | its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch |
|---|
| 149 | jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch |
|---|
| 150 | CASE ( DATA_ORDER_XYZ ) |
|---|
| 151 | ids = grid%sd31 ; ide = grid%ed31 ; |
|---|
| 152 | jds = grid%sd32 ; jde = grid%ed32 ; |
|---|
| 153 | kds = grid%sd33 ; kde = grid%ed33 ; |
|---|
| 154 | |
|---|
| 155 | ims = grid%sm31 ; ime = grid%em31 ; |
|---|
| 156 | jms = grid%sm32 ; jme = grid%em32 ; |
|---|
| 157 | kms = grid%sm33 ; kme = grid%em33 ; |
|---|
| 158 | |
|---|
| 159 | its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch |
|---|
| 160 | jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch |
|---|
| 161 | kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch |
|---|
| 162 | CASE ( DATA_ORDER_XZY ) |
|---|
| 163 | ids = grid%sd31 ; ide = grid%ed31 ; |
|---|
| 164 | kds = grid%sd32 ; kde = grid%ed32 ; |
|---|
| 165 | jds = grid%sd33 ; jde = grid%ed33 ; |
|---|
| 166 | |
|---|
| 167 | ims = grid%sm31 ; ime = grid%em31 ; |
|---|
| 168 | kms = grid%sm32 ; kme = grid%em32 ; |
|---|
| 169 | jms = grid%sm33 ; jme = grid%em33 ; |
|---|
| 170 | |
|---|
| 171 | its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch |
|---|
| 172 | kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch |
|---|
| 173 | jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch |
|---|
| 174 | |
|---|
| 175 | END SELECT |
|---|
| 176 | |
|---|
| 177 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 178 | !!MARS : mountain |
|---|
| 179 | !!MARS : mountain ex. hm = 2000. xa = 6.0 |
|---|
| 180 | open(unit=22,file='ze_hill',form='formatted',status='old') |
|---|
| 181 | rewind(22) |
|---|
| 182 | read(22,*) hm, xa |
|---|
| 183 | write(6,*) 'height, width ', hm, xa |
|---|
| 184 | close(22) |
|---|
| 185 | !!MARS |
|---|
| 186 | !!MARS |
|---|
| 187 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 188 | |
|---|
| 189 | |
|---|
| 190 | icm = ide/2 |
|---|
| 191 | !****Mars |
|---|
| 192 | jcm = jde/2 |
|---|
| 193 | ! |
|---|
| 194 | ! xa1 = 5000./500. |
|---|
| 195 | ! xal1 = 4000./500. |
|---|
| 196 | ! pii = 2.*asin(1.0) |
|---|
| 197 | ! hm1 = 250. |
|---|
| 198 | !! hm1 = 1000. |
|---|
| 199 | !****Mars |
|---|
| 200 | |
|---|
| 201 | delt = 3. |
|---|
| 202 | ! delt = 10. |
|---|
| 203 | |
|---|
| 204 | !****Mars |
|---|
| 205 | stretch_grid = .true. |
|---|
| 206 | ! stretch_grid = .false. |
|---|
| 207 | !****Mars |
|---|
| 208 | ! z_scale = .50 |
|---|
| 209 | z_scale = .40 |
|---|
| 210 | pi = 2.*asin(1.0) |
|---|
| 211 | write(6,*) ' pi is ',pi |
|---|
| 212 | nxc = (ide-ids)/2 |
|---|
| 213 | nyc = (jde-jds)/2 |
|---|
| 214 | |
|---|
| 215 | !!!MARS |
|---|
| 216 | !!!MARS |
|---|
| 217 | ! open(unit=16,file='input_vert',form='formatted',status='old') |
|---|
| 218 | ! rewind(16) |
|---|
| 219 | ! read(16,*) delt, z_scale |
|---|
| 220 | ! write(6,*) 'delt, z_scale are ', delt, z_scale |
|---|
| 221 | ! close(16) |
|---|
| 222 | !!!MARS |
|---|
| 223 | !!!MARS |
|---|
| 224 | |
|---|
| 225 | CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) |
|---|
| 226 | |
|---|
| 227 | ! here we check to see if the boundary conditions are set properly |
|---|
| 228 | |
|---|
| 229 | CALL boundary_condition_check( config_flags, bdyzone, error, grid%id ) |
|---|
| 230 | |
|---|
| 231 | moisture_init = .true. |
|---|
| 232 | |
|---|
| 233 | grid%itimestep=0 |
|---|
| 234 | |
|---|
| 235 | #ifdef DM_PARALLEL |
|---|
| 236 | CALL wrf_dm_bcast_bytes( icm , IWORDSIZE ) |
|---|
| 237 | CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE ) |
|---|
| 238 | #endif |
|---|
| 239 | |
|---|
| 240 | CALL nl_set_mminlu(1,' ') |
|---|
| 241 | CALL nl_set_iswater(1,0) |
|---|
| 242 | CALL nl_set_cen_lat(1,40.) |
|---|
| 243 | CALL nl_set_cen_lon(1,-105.) |
|---|
| 244 | CALL nl_set_truelat1(1,0.) |
|---|
| 245 | CALL nl_set_truelat2(1,0.) |
|---|
| 246 | CALL nl_set_moad_cen_lat (1,0.) |
|---|
| 247 | CALL nl_set_stand_lon (1,0.) |
|---|
| 248 | CALL nl_set_map_proj(1,0) |
|---|
| 249 | |
|---|
| 250 | |
|---|
| 251 | ! here we initialize data we currently is not initialized |
|---|
| 252 | ! in the input data |
|---|
| 253 | |
|---|
| 254 | DO j = jts, jte |
|---|
| 255 | DO i = its, ite |
|---|
| 256 | grid%msft(i,j) = 1. |
|---|
| 257 | grid%msfu(i,j) = 1. |
|---|
| 258 | grid%msfv(i,j) = 1. |
|---|
| 259 | grid%sina(i,j) = 0. |
|---|
| 260 | grid%cosa(i,j) = 1. |
|---|
| 261 | grid%e(i,j) = 0. |
|---|
| 262 | grid%f(i,j) = 0. !! MARS: put coriolis here if needed |
|---|
| 263 | |
|---|
| 264 | END DO |
|---|
| 265 | END DO |
|---|
| 266 | |
|---|
| 267 | DO j = jts, jte |
|---|
| 268 | DO k = kts, kte |
|---|
| 269 | DO i = its, ite |
|---|
| 270 | grid%em_ww(i,k,j) = 0. |
|---|
| 271 | END DO |
|---|
| 272 | END DO |
|---|
| 273 | END DO |
|---|
| 274 | |
|---|
| 275 | grid%step_number = 0 |
|---|
| 276 | |
|---|
| 277 | !! set up the grid |
|---|
| 278 | ! |
|---|
| 279 | ! IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz) |
|---|
| 280 | ! DO k=1, kde |
|---|
| 281 | ! grid%em_znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ & |
|---|
| 282 | ! (1.-exp(-1./z_scale)) |
|---|
| 283 | ! ENDDO |
|---|
| 284 | ! ELSE |
|---|
| 285 | ! DO k=1, kde |
|---|
| 286 | ! grid%em_znw(k) = 1. - float(k-1)/float(kde-1) |
|---|
| 287 | ! ENDDO |
|---|
| 288 | ! ENDIF |
|---|
| 289 | |
|---|
| 290 | !!MARS |
|---|
| 291 | !!MARS |
|---|
| 292 | open(unit=12,file='levels',form='formatted',status='old') |
|---|
| 293 | rewind(12) |
|---|
| 294 | DO k=1, kde |
|---|
| 295 | read(12,*) grid%em_znw(k) |
|---|
| 296 | write(6,*) 'read level ', k,grid%em_znw(k) |
|---|
| 297 | ENDDO |
|---|
| 298 | close(12) |
|---|
| 299 | !!MARS |
|---|
| 300 | !!MARS |
|---|
| 301 | |
|---|
| 302 | |
|---|
| 303 | DO k=1, kde-1 |
|---|
| 304 | grid%em_dnw(k) = grid%em_znw(k+1) - grid%em_znw(k) |
|---|
| 305 | grid%em_rdnw(k) = 1./grid%em_dnw(k) |
|---|
| 306 | grid%em_znu(k) = 0.5*(grid%em_znw(k+1)+grid%em_znw(k)) |
|---|
| 307 | ENDDO |
|---|
| 308 | DO k=2, kde-1 |
|---|
| 309 | grid%em_dn(k) = 0.5*(grid%em_dnw(k)+grid%em_dnw(k-1)) |
|---|
| 310 | grid%em_rdn(k) = 1./grid%em_dn(k) |
|---|
| 311 | grid%em_fnp(k) = .5* grid%em_dnw(k )/grid%em_dn(k) |
|---|
| 312 | grid%em_fnm(k) = .5* grid%em_dnw(k-1)/grid%em_dn(k) |
|---|
| 313 | ENDDO |
|---|
| 314 | |
|---|
| 315 | cof1 = (2.*grid%em_dn(2)+grid%em_dn(3))/(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(2) |
|---|
| 316 | cof2 = grid%em_dn(2) /(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(3) |
|---|
| 317 | grid%cf1 = grid%em_fnp(2) + cof1 |
|---|
| 318 | grid%cf2 = grid%em_fnm(2) - cof1 - cof2 |
|---|
| 319 | grid%cf3 = cof2 |
|---|
| 320 | |
|---|
| 321 | grid%cfn = (.5*grid%em_dnw(kde-1)+grid%em_dn(kde-1))/grid%em_dn(kde-1) |
|---|
| 322 | grid%cfn1 = -.5*grid%em_dnw(kde-1)/grid%em_dn(kde-1) |
|---|
| 323 | grid%rdx = 1./config_flags%dx |
|---|
| 324 | grid%rdy = 1./config_flags%dy |
|---|
| 325 | |
|---|
| 326 | ! get the sounding from the ascii sounding file, first get dry sounding and |
|---|
| 327 | ! calculate base state |
|---|
| 328 | |
|---|
| 329 | !!!! |
|---|
| 330 | !!!! user-modified wind speed |
|---|
| 331 | !!!! |
|---|
| 332 | mulu = 1. !! default |
|---|
| 333 | mulv = 1. !! default |
|---|
| 334 | addu = 0. !! default |
|---|
| 335 | addv = 0. !! default |
|---|
| 336 | IF (config_flags%init_MU .ne. 0.) mulu = config_flags%init_MU |
|---|
| 337 | IF (config_flags%init_MV .ne. 0.) mulv = config_flags%init_MV |
|---|
| 338 | IF (config_flags%init_U .ne. 0.) addu = config_flags%init_U |
|---|
| 339 | IF (config_flags%init_V .ne. 0.) addv = config_flags%init_V |
|---|
| 340 | write(6,*) ' coeff for winds: ', mulu, mulv, addu, addv |
|---|
| 341 | |
|---|
| 342 | dry_sounding = .true. |
|---|
| 343 | IF ( wrf_dm_on_monitor() ) THEN |
|---|
| 344 | write(6,*) ' getting dry sounding for base state ' |
|---|
| 345 | |
|---|
| 346 | CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in, & |
|---|
| 347 | mulu, mulv, addu, addv ) |
|---|
| 348 | ENDIF |
|---|
| 349 | CALL wrf_dm_bcast_real( zk , nl_max ) |
|---|
| 350 | CALL wrf_dm_bcast_real( p_in , nl_max ) |
|---|
| 351 | CALL wrf_dm_bcast_real( pd_in , nl_max ) |
|---|
| 352 | CALL wrf_dm_bcast_real( theta , nl_max ) |
|---|
| 353 | CALL wrf_dm_bcast_real( rho , nl_max ) |
|---|
| 354 | CALL wrf_dm_bcast_real( u , nl_max ) |
|---|
| 355 | CALL wrf_dm_bcast_real( v , nl_max ) |
|---|
| 356 | CALL wrf_dm_bcast_real( qv , nl_max ) |
|---|
| 357 | CALL wrf_dm_bcast_integer ( nl_in , 1 ) |
|---|
| 358 | |
|---|
| 359 | write(6,*) ' returned from reading sounding, nl_in is ',nl_in |
|---|
| 360 | |
|---|
| 361 | ! find ptop for the desired ztop (ztop is input from the namelist), |
|---|
| 362 | ! and find surface pressure |
|---|
| 363 | |
|---|
| 364 | grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in ) |
|---|
| 365 | |
|---|
| 366 | !!MARS |
|---|
| 367 | !!MARS |
|---|
| 368 | open(unit=14,file='input_coord',form='formatted',status='old') |
|---|
| 369 | rewind(14) |
|---|
| 370 | read(14,*) lon_input |
|---|
| 371 | read(14,*) lat_input |
|---|
| 372 | close(14) |
|---|
| 373 | write(6,*) ' lon is ',lon_input |
|---|
| 374 | write(6,*) ' lat is ',lat_input |
|---|
| 375 | !!MARS |
|---|
| 376 | !!MARS |
|---|
| 377 | |
|---|
| 378 | !!MARS |
|---|
| 379 | !!MARS |
|---|
| 380 | open(unit=18,file='input_more',form='formatted',status='old') |
|---|
| 381 | rewind(18) |
|---|
| 382 | read(18,*) alt_input, tsurf_input |
|---|
| 383 | close(18) |
|---|
| 384 | write(6,*) ' alt is ',alt_input |
|---|
| 385 | write(6,*) ' tsurf is ',tsurf_input |
|---|
| 386 | !!MARS |
|---|
| 387 | !!MARS |
|---|
| 388 | |
|---|
| 389 | DO j=jts,jte |
|---|
| 390 | DO i=its,ite |
|---|
| 391 | !!MARS |
|---|
| 392 | grid%ht(i,j) = alt_input |
|---|
| 393 | ! flat surface |
|---|
| 394 | !! grid%ht(i,j) = 0. |
|---|
| 395 | ! grid%ht(i,j) = hm1*exp(-(( float(i-icm)/xa1)**2)) & |
|---|
| 396 | ! *( (cos(pii*float(i-icm)/xal1))**2 ) |
|---|
| 397 | !****Mars |
|---|
| 398 | !!3D hill |
|---|
| 399 | ! grid%ht(i,j) = hm/(1.+(float(i-icm)/xa)**2+(float(j-jcm)/xa)**2) |
|---|
| 400 | |
|---|
| 401 | !!2D hill |
|---|
| 402 | grid%ht(i,j) = alt_input + hm/(1.+(float(i-icm)/xa)**2) |
|---|
| 403 | |
|---|
| 404 | !!!3D crater |
|---|
| 405 | !! grid%ht(i,j) = hm - hm/(1.+(float(i-icm)/xa)**2+(float(j-jcm)/xa)**2) |
|---|
| 406 | !!3D crater w/ rims |
|---|
| 407 | ! x_param = float(i-icm) |
|---|
| 408 | ! y_param = float(j-jcm) |
|---|
| 409 | ! dilat = xa/2 |
|---|
| 410 | ! rho_param = sqrt(x_param**2 + y_param**2) |
|---|
| 411 | ! ! revolution surface ; seed is a fourth order polynom |
|---|
| 412 | ! grid%ht(i,j) = (rho_param+6*dilat)*(rho_param+10*dilat) |
|---|
| 413 | ! grid%ht(i,j) = (rho_param-6*dilat)*(rho_param-10*dilat)*grid%ht(i,j) |
|---|
| 414 | ! ! flat terrain elsewhere - smooth gradient (no abrupt fall) |
|---|
| 415 | ! grid%ht(i,j) = grid%ht(i,j)*(tanh(rho_param+7*dilat)/2 - tanh(rho_param-7*dilat)/2) |
|---|
| 416 | ! grid%ht(i,j) = hm - (hm*.4/1500)*grid%ht(i,j)/(dilat**4) |
|---|
| 417 | ! !NONONONONON |
|---|
| 418 | ! !grid%ht(i,j) = grid%ht(i,j) + alt_input |
|---|
| 419 | ! !if (rho_param .GE. dilat*10) ht(i,j) = hm |
|---|
| 420 | |
|---|
| 421 | grid%tsk(i,j) = tsurf_input |
|---|
| 422 | !!MARS |
|---|
| 423 | grid%xlat(i,j) = lat_input |
|---|
| 424 | grid%xlong(i,j) = lon_input!+float(i)*config_flags%dx/59000. |
|---|
| 425 | grid%mars_emiss(i,j)=0.95 |
|---|
| 426 | grid%mars_cice(i,j)=0. |
|---|
| 427 | grid%slpx(i,j) = 0. |
|---|
| 428 | grid%slpy(i,j) = 0. |
|---|
| 429 | !!MARS |
|---|
| 430 | ENDDO |
|---|
| 431 | ENDDO |
|---|
| 432 | |
|---|
| 433 | xs=ide/2 -3 |
|---|
| 434 | xs=ids -3 |
|---|
| 435 | xe=xs + 6 |
|---|
| 436 | ys=jde/2 -3 |
|---|
| 437 | ye=ys + 6 |
|---|
| 438 | mtn_ht = 500 |
|---|
| 439 | #ifdef MTN |
|---|
| 440 | DO j=max(ys,jds),min(ye,jde-1) |
|---|
| 441 | DO i=max(xs,ids),min(xe,ide-1) |
|---|
| 442 | grid%ht(i,j) = mtn_ht * 0.25 * & |
|---|
| 443 | ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) * & |
|---|
| 444 | ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) ) |
|---|
| 445 | ENDDO |
|---|
| 446 | ENDDO |
|---|
| 447 | #endif |
|---|
| 448 | #ifdef EW_RIDGE |
|---|
| 449 | DO j=max(ys,jds),min(ye,jde-1) |
|---|
| 450 | DO i=ids,ide |
|---|
| 451 | grid%ht(i,j) = mtn_ht * 0.50 * & |
|---|
| 452 | ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) ) |
|---|
| 453 | ENDDO |
|---|
| 454 | ENDDO |
|---|
| 455 | #endif |
|---|
| 456 | #ifdef NS_RIDGE |
|---|
| 457 | DO j=jds,jde |
|---|
| 458 | DO i=max(xs,ids),min(xe,ide-1) |
|---|
| 459 | grid%ht(i,j) = mtn_ht * 0.50 * & |
|---|
| 460 | ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) |
|---|
| 461 | ENDDO |
|---|
| 462 | ENDDO |
|---|
| 463 | #endif |
|---|
| 464 | DO j=jts,jte |
|---|
| 465 | DO i=its,ite |
|---|
| 466 | grid%em_phb(i,1,j) = g * grid%ht(i,j) |
|---|
| 467 | grid%em_ph0(i,1,j) = g * grid%ht(i,j) |
|---|
| 468 | ENDDO |
|---|
| 469 | ENDDO |
|---|
| 470 | |
|---|
| 471 | !!!dans hill_2d OK |
|---|
| 472 | ! grid%em_phb(i,1,j) = g*grid%ht(i,j) |
|---|
| 473 | ! grid%em_php(i,1,j) = 0. |
|---|
| 474 | ! grid%em_ph0(i,1,j) = grid%em_phb(i,1,j) |
|---|
| 475 | |
|---|
| 476 | DO J = jts, jte |
|---|
| 477 | DO I = its, ite |
|---|
| 478 | |
|---|
| 479 | p_surf = interp_0( p_in, zk, grid%em_phb(i,1,j)/g, nl_in ) |
|---|
| 480 | grid%em_mub(i,j) = p_surf-grid%p_top |
|---|
| 481 | |
|---|
| 482 | ! this is dry hydrostatic sounding (base state), so given grid%em_p (coordinate), |
|---|
| 483 | ! interp theta (from interp) and compute 1/rho from eqn. of state |
|---|
| 484 | |
|---|
| 485 | DO K = 1, kte-1 |
|---|
| 486 | p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top |
|---|
| 487 | grid%em_pb(i,k,j) = p_level |
|---|
| 488 | grid%em_t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0 |
|---|
| 489 | grid%em_alb(i,k,j) = (r_d/p1000mb)*(grid%em_t_init(i,k,j)+t0)*(grid%em_pb(i,k,j)/p1000mb)**cvpm |
|---|
| 490 | ENDDO |
|---|
| 491 | |
|---|
| 492 | ! calc hydrostatic balance (alternatively we could interp the geopotential from the |
|---|
| 493 | ! sounding, but this assures that the base state is in exact hydrostatic balance with |
|---|
| 494 | ! respect to the model eqns. |
|---|
| 495 | |
|---|
| 496 | DO k = 2,kte |
|---|
| 497 | grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j) |
|---|
| 498 | ENDDO |
|---|
| 499 | |
|---|
| 500 | ENDDO |
|---|
| 501 | ENDDO |
|---|
| 502 | |
|---|
| 503 | IF ( wrf_dm_on_monitor() ) THEN |
|---|
| 504 | write(6,*) ' ptop is ',grid%p_top |
|---|
| 505 | write(6,*) ' base state grid%em_mub(1,1), p_surf is ',grid%em_mub(1,1),grid%em_mub(1,1)+grid%p_top |
|---|
| 506 | ENDIF |
|---|
| 507 | |
|---|
| 508 | ! calculate full state for each column - this includes moisture. |
|---|
| 509 | |
|---|
| 510 | !!!!!MARS MARS |
|---|
| 511 | ! write(6,*) ' getting moist sounding for full state ' |
|---|
| 512 | ! dry_sounding = .false. |
|---|
| 513 | dry_sounding = .true. |
|---|
| 514 | CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in, & |
|---|
| 515 | mulu, mulv, addu, addv ) |
|---|
| 516 | |
|---|
| 517 | DO J = jts, min(jde-1,jte) |
|---|
| 518 | DO I = its, min(ide-1,ite) |
|---|
| 519 | |
|---|
| 520 | ! At this point grid%p_top is already set. find the DRY mass in the column |
|---|
| 521 | ! by interpolating the DRY pressure. |
|---|
| 522 | |
|---|
| 523 | pd_surf = interp_0( pd_in, zk, grid%em_phb(i,1,j)/g, nl_in ) |
|---|
| 524 | |
|---|
| 525 | ! compute the perturbation mass and the full mass |
|---|
| 526 | |
|---|
| 527 | grid%em_mu_1(i,j) = pd_surf-grid%p_top - grid%em_mub(i,j) |
|---|
| 528 | grid%em_mu_2(i,j) = grid%em_mu_1(i,j) |
|---|
| 529 | grid%em_mu0(i,j) = grid%em_mu_1(i,j) + grid%em_mub(i,j) |
|---|
| 530 | |
|---|
| 531 | ! given the dry pressure and coordinate system, interp the potential |
|---|
| 532 | ! temperature and qv |
|---|
| 533 | |
|---|
| 534 | do k=1,kde-1 |
|---|
| 535 | |
|---|
| 536 | p_level = grid%em_znu(k)*(pd_surf - grid%p_top) + grid%p_top |
|---|
| 537 | |
|---|
| 538 | moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in ) |
|---|
| 539 | grid%em_t_1(i,k,j) = interp_0( theta, pd_in, p_level, nl_in ) - t0 |
|---|
| 540 | grid%em_t_2(i,k,j) = grid%em_t_1(i,k,j) |
|---|
| 541 | |
|---|
| 542 | |
|---|
| 543 | enddo |
|---|
| 544 | |
|---|
| 545 | ! integrate the hydrostatic equation (from the RHS of the bigstep |
|---|
| 546 | ! vertical momentum equation) down from the top to get grid%em_p. |
|---|
| 547 | ! first from the top of the model to the top pressure |
|---|
| 548 | |
|---|
| 549 | k = kte-1 ! top level |
|---|
| 550 | |
|---|
| 551 | qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV)) |
|---|
| 552 | qvf2 = 1./(1.+qvf1) |
|---|
| 553 | qvf1 = qvf1*qvf2 |
|---|
| 554 | |
|---|
| 555 | ! grid%em_p(i,k,j) = - 0.5*grid%em_mu_1(i,j)/grid%em_rdnw(k) |
|---|
| 556 | grid%em_p(i,k,j) = - 0.5*(grid%em_mu_1(i,j)+qvf1*grid%em_mub(i,j))/grid%em_rdnw(k)/qvf2 |
|---|
| 557 | qvf = 1. + rvovrd*moist(i,k,j,P_QV) |
|---|
| 558 | grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* & |
|---|
| 559 | (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm) |
|---|
| 560 | grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j) |
|---|
| 561 | |
|---|
| 562 | ! down the column |
|---|
| 563 | |
|---|
| 564 | do k=kte-2,1,-1 |
|---|
| 565 | qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV)) |
|---|
| 566 | qvf2 = 1./(1.+qvf1) |
|---|
| 567 | qvf1 = qvf1*qvf2 |
|---|
| 568 | grid%em_p(i,k,j) = grid%em_p(i,k+1,j) - (grid%em_mu_1(i,j) + qvf1*grid%em_mub(i,j))/qvf2/grid%em_rdn(k+1) |
|---|
| 569 | qvf = 1. + rvovrd*moist(i,k,j,P_QV) |
|---|
| 570 | grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* & |
|---|
| 571 | (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm) |
|---|
| 572 | grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j) |
|---|
| 573 | enddo |
|---|
| 574 | |
|---|
| 575 | ! this is the hydrostatic equation used in the model after the |
|---|
| 576 | ! small timesteps. In the model, grid%em_al (inverse density) |
|---|
| 577 | ! is computed from the geopotential. |
|---|
| 578 | |
|---|
| 579 | |
|---|
| 580 | grid%em_ph_1(i,1,j) = 0. |
|---|
| 581 | DO k = 2,kte |
|---|
| 582 | grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*( & |
|---|
| 583 | (grid%em_mub(i,j)+grid%em_mu_1(i,j))*grid%em_al(i,k-1,j)+ & |
|---|
| 584 | grid%em_mu_1(i,j)*grid%em_alb(i,k-1,j) ) |
|---|
| 585 | |
|---|
| 586 | grid%em_ph_2(i,k,j) = grid%em_ph_1(i,k,j) |
|---|
| 587 | grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j) |
|---|
| 588 | ENDDO |
|---|
| 589 | |
|---|
| 590 | IF ( wrf_dm_on_monitor() ) THEN |
|---|
| 591 | if((i==2) .and. (j==2)) then |
|---|
| 592 | write(6,*) ' grid%em_ph_1 calc ',grid%em_ph_1(2,1,2),grid%em_ph_1(2,2,2),& |
|---|
| 593 | grid%em_mu_1(2,2)+grid%em_mub(2,2),grid%em_mu_1(2,2), & |
|---|
| 594 | grid%em_alb(2,1,2),grid%em_al(1,2,1),grid%em_rdnw(1) |
|---|
| 595 | endif |
|---|
| 596 | ENDIF |
|---|
| 597 | |
|---|
| 598 | ENDDO |
|---|
| 599 | ENDDO |
|---|
| 600 | |
|---|
| 601 | !#if 0 |
|---|
| 602 | |
|---|
| 603 | ! thermal perturbation to kick off convection |
|---|
| 604 | |
|---|
| 605 | write(6,*) ' nxc, nyc for perturbation ',nxc,nyc |
|---|
| 606 | write(6,*) ' delt for perturbation ',delt |
|---|
| 607 | |
|---|
| 608 | DO J = jts, min(jde-1,jte) |
|---|
| 609 | yrad = config_flags%dy*float(j-nyc)/10000. |
|---|
| 610 | yrad = 0. |
|---|
| 611 | DO I = its, min(ide-1,ite) |
|---|
| 612 | xrad = config_flags%dx*float(i-nxc)/10000. |
|---|
| 613 | xrad = 0. |
|---|
| 614 | DO K = 1, kte-1 |
|---|
| 615 | |
|---|
| 616 | ! put in preturbation theta (bubble) and recalc density. note, |
|---|
| 617 | ! the mass in the column is not changing, so when theta changes, |
|---|
| 618 | ! we recompute density and geopotential |
|---|
| 619 | |
|---|
| 620 | zrad = 0.5*(grid%em_ph_1(i,k,j)+grid%em_ph_1(i,k+1,j) & |
|---|
| 621 | +grid%em_phb(i,k,j)+grid%em_phb(i,k+1,j))/g |
|---|
| 622 | zrad = (zrad-1500.)/1500. |
|---|
| 623 | RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad) |
|---|
| 624 | IF(RAD <= 1.) THEN |
|---|
| 625 | grid%em_t_1(i,k,j)=grid%em_t_1(i,k,j)+delt*COS(.5*PI*RAD)**2 |
|---|
| 626 | grid%em_t_2(i,k,j)=grid%em_t_1(i,k,j) |
|---|
| 627 | qvf = 1. + rvovrd*moist(i,k,j,P_QV) |
|---|
| 628 | grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* & |
|---|
| 629 | (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm) |
|---|
| 630 | grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j) |
|---|
| 631 | ENDIF |
|---|
| 632 | ENDDO |
|---|
| 633 | |
|---|
| 634 | ! rebalance hydrostatically |
|---|
| 635 | |
|---|
| 636 | DO k = 2,kte |
|---|
| 637 | grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*( & |
|---|
| 638 | (grid%em_mub(i,j)+grid%em_mu_1(i,j))*grid%em_al(i,k-1,j)+ & |
|---|
| 639 | grid%em_mu_1(i,j)*grid%em_alb(i,k-1,j) ) |
|---|
| 640 | |
|---|
| 641 | grid%em_ph_2(i,k,j) = grid%em_ph_1(i,k,j) |
|---|
| 642 | grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j) |
|---|
| 643 | ENDDO |
|---|
| 644 | |
|---|
| 645 | ENDDO |
|---|
| 646 | ENDDO |
|---|
| 647 | |
|---|
| 648 | !#endif |
|---|
| 649 | |
|---|
| 650 | IF ( wrf_dm_on_monitor() ) THEN |
|---|
| 651 | write(6,*) ' grid%em_mu_1 from comp ', grid%em_mu_1(1,1) |
|---|
| 652 | write(6,*) ' full state sounding from comp, ph, grid%em_p, grid%em_al, grid%em_t_1, qv ' |
|---|
| 653 | do k=1,kde-1 |
|---|
| 654 | write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%em_ph_1(1,k,1)+grid%em_phb(1,k,1), & |
|---|
| 655 | grid%em_p(1,k,1)+grid%em_pb(1,k,1), grid%em_alt(1,k,1), & |
|---|
| 656 | grid%em_t_1(1,k,1)+t0, moist(1,k,1,P_QV) |
|---|
| 657 | enddo |
|---|
| 658 | |
|---|
| 659 | write(6,*) ' pert state sounding from comp, grid%em_ph_1, pp, alp, grid%em_t_1, qv ' |
|---|
| 660 | do k=1,kde-1 |
|---|
| 661 | write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%em_ph_1(1,k,1), & |
|---|
| 662 | grid%em_p(1,k,1), grid%em_al(1,k,1), & |
|---|
| 663 | grid%em_t_1(1,k,1), moist(1,k,1,P_QV) |
|---|
| 664 | enddo |
|---|
| 665 | ENDIF |
|---|
| 666 | |
|---|
| 667 | ! interp v |
|---|
| 668 | |
|---|
| 669 | DO J = jts, jte |
|---|
| 670 | DO I = its, min(ide-1,ite) |
|---|
| 671 | |
|---|
| 672 | IF (j == jds) THEN |
|---|
| 673 | z_at_v = grid%em_phb(i,1,j)/g |
|---|
| 674 | ELSE IF (j == jde) THEN |
|---|
| 675 | z_at_v = grid%em_phb(i,1,j-1)/g |
|---|
| 676 | ELSE |
|---|
| 677 | z_at_v = 0.5*(grid%em_phb(i,1,j)+grid%em_phb(i,1,j-1))/g |
|---|
| 678 | END IF |
|---|
| 679 | |
|---|
| 680 | p_surf = interp_0( p_in, zk, z_at_v, nl_in ) |
|---|
| 681 | |
|---|
| 682 | DO K = 1, kte-1 |
|---|
| 683 | p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top |
|---|
| 684 | grid%em_v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in ) |
|---|
| 685 | grid%em_v_2(i,k,j) = grid%em_v_1(i,k,j) |
|---|
| 686 | ENDDO |
|---|
| 687 | |
|---|
| 688 | ENDDO |
|---|
| 689 | ENDDO |
|---|
| 690 | |
|---|
| 691 | ! interp u |
|---|
| 692 | |
|---|
| 693 | DO J = jts, min(jde-1,jte) |
|---|
| 694 | DO I = its, ite |
|---|
| 695 | |
|---|
| 696 | IF (i == ids) THEN |
|---|
| 697 | z_at_u = grid%em_phb(i,1,j)/g |
|---|
| 698 | ELSE IF (i == ide) THEN |
|---|
| 699 | z_at_u = grid%em_phb(i-1,1,j)/g |
|---|
| 700 | ELSE |
|---|
| 701 | z_at_u = 0.5*(grid%em_phb(i,1,j)+grid%em_phb(i-1,1,j))/g |
|---|
| 702 | END IF |
|---|
| 703 | |
|---|
| 704 | p_surf = interp_0( p_in, zk, z_at_u, nl_in ) |
|---|
| 705 | |
|---|
| 706 | DO K = 1, kte-1 |
|---|
| 707 | p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top |
|---|
| 708 | grid%em_u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in ) |
|---|
| 709 | grid%em_u_2(i,k,j) = grid%em_u_1(i,k,j) |
|---|
| 710 | ENDDO |
|---|
| 711 | |
|---|
| 712 | ENDDO |
|---|
| 713 | ENDDO |
|---|
| 714 | |
|---|
| 715 | ! set w |
|---|
| 716 | |
|---|
| 717 | DO J = jts, min(jde-1,jte) |
|---|
| 718 | DO K = kts, kte |
|---|
| 719 | DO I = its, min(ide-1,ite) |
|---|
| 720 | grid%em_w_1(i,k,j) = 0. |
|---|
| 721 | grid%em_w_2(i,k,j) = 0. |
|---|
| 722 | ENDDO |
|---|
| 723 | ENDDO |
|---|
| 724 | ENDDO |
|---|
| 725 | |
|---|
| 726 | ! set a few more things |
|---|
| 727 | |
|---|
| 728 | DO J = jts, min(jde-1,jte) |
|---|
| 729 | DO K = kts, kte-1 |
|---|
| 730 | DO I = its, min(ide-1,ite) |
|---|
| 731 | grid%h_diabatic(i,k,j) = 0. |
|---|
| 732 | ENDDO |
|---|
| 733 | ENDDO |
|---|
| 734 | ENDDO |
|---|
| 735 | |
|---|
| 736 | IF ( wrf_dm_on_monitor() ) THEN |
|---|
| 737 | DO k=1,kte-1 |
|---|
| 738 | grid%em_t_base(k) = grid%em_t_1(1,k,1) |
|---|
| 739 | grid%qv_base(k) = moist(1,k,1,P_QV) |
|---|
| 740 | grid%u_base(k) = grid%em_u_1(1,k,1) |
|---|
| 741 | grid%v_base(k) = grid%em_v_1(1,k,1) |
|---|
| 742 | grid%z_base(k) = 0.5*(grid%em_phb(1,k,1)+grid%em_phb(1,k+1,1)+grid%em_ph_1(1,k,1)+grid%em_ph_1(1,k+1,1))/g |
|---|
| 743 | ENDDO |
|---|
| 744 | ENDIF |
|---|
| 745 | CALL wrf_dm_bcast_real( grid%em_t_base , kte ) |
|---|
| 746 | CALL wrf_dm_bcast_real( grid%qv_base , kte ) |
|---|
| 747 | CALL wrf_dm_bcast_real( grid%u_base , kte ) |
|---|
| 748 | CALL wrf_dm_bcast_real( grid%v_base , kte ) |
|---|
| 749 | CALL wrf_dm_bcast_real( grid%z_base , kte ) |
|---|
| 750 | |
|---|
| 751 | DO J = jts, min(jde-1,jte) |
|---|
| 752 | DO I = its, min(ide-1,ite) |
|---|
| 753 | thtmp = grid%em_t_2(i,1,j)+t0 |
|---|
| 754 | ptmp = grid%em_p(i,1,j)+grid%em_pb(i,1,j) |
|---|
| 755 | temp(1) = thtmp * (ptmp/p1000mb)**rcp |
|---|
| 756 | thtmp = grid%em_t_2(i,2,j)+t0 |
|---|
| 757 | ptmp = grid%em_p(i,2,j)+grid%em_pb(i,2,j) |
|---|
| 758 | temp(2) = thtmp * (ptmp/p1000mb)**rcp |
|---|
| 759 | thtmp = grid%em_t_2(i,3,j)+t0 |
|---|
| 760 | ptmp = grid%em_p(i,3,j)+grid%em_pb(i,3,j) |
|---|
| 761 | temp(3) = thtmp * (ptmp/p1000mb)**rcp |
|---|
| 762 | |
|---|
| 763 | !!MARS |
|---|
| 764 | ! grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3) |
|---|
| 765 | grid%tmn(I,J)=grid%tsk(I,J)-0.5 |
|---|
| 766 | !!!MARS |
|---|
| 767 | !!TODO: passer la valeur a partir des donnees |
|---|
| 768 | !grid%mars_tsoil(I,:,J)=grid%tsk(I,J) |
|---|
| 769 | !!!MARS |
|---|
| 770 | ENDDO |
|---|
| 771 | ENDDO |
|---|
| 772 | |
|---|
| 773 | END SUBROUTINE init_domain_rk |
|---|
| 774 | |
|---|
| 775 | SUBROUTINE init_module_initialize |
|---|
| 776 | END SUBROUTINE init_module_initialize |
|---|
| 777 | |
|---|
| 778 | !--------------------------------------------------------------------- |
|---|
| 779 | |
|---|
| 780 | ! test driver for get_sounding |
|---|
| 781 | ! |
|---|
| 782 | ! implicit none |
|---|
| 783 | ! integer n |
|---|
| 784 | ! parameter(n = 1000) |
|---|
| 785 | ! real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n) |
|---|
| 786 | ! logical dry |
|---|
| 787 | ! integer nl,k |
|---|
| 788 | ! |
|---|
| 789 | ! dry = .false. |
|---|
| 790 | ! dry = .true. |
|---|
| 791 | ! call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl ) |
|---|
| 792 | ! write(6,*) ' input levels ',nl |
|---|
| 793 | ! write(6,*) ' sounding ' |
|---|
| 794 | ! write(6,*) ' k height(m) press (Pa) pd(Pa) theta (K) den(kg/m^3) u(m/s) v(m/s) qv(g/g) ' |
|---|
| 795 | ! do k=1,nl |
|---|
| 796 | ! 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) |
|---|
| 797 | ! enddo |
|---|
| 798 | ! end |
|---|
| 799 | ! |
|---|
| 800 | !--------------------------------------------------------------------------- |
|---|
| 801 | |
|---|
| 802 | subroutine get_sounding( zk, p, p_dry, theta, rho, & |
|---|
| 803 | u, v, qv, dry, nl_max, nl_in, & |
|---|
| 804 | mulu, mulv, addu, addv ) |
|---|
| 805 | implicit none |
|---|
| 806 | |
|---|
| 807 | integer nl_max, nl_in |
|---|
| 808 | real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), & |
|---|
| 809 | u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max) |
|---|
| 810 | logical dry |
|---|
| 811 | |
|---|
| 812 | integer n |
|---|
| 813 | parameter(n=1000) |
|---|
| 814 | logical debug |
|---|
| 815 | |
|---|
| 816 | ! parameter( debug = .false.) |
|---|
| 817 | !****Mars |
|---|
| 818 | parameter( debug = .true.) |
|---|
| 819 | real mulu, mulv, addu, addv |
|---|
| 820 | |
|---|
| 821 | |
|---|
| 822 | ! input sounding data |
|---|
| 823 | |
|---|
| 824 | real p_surf, th_surf, qv_surf |
|---|
| 825 | real pi_surf, pi(n) |
|---|
| 826 | real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n) |
|---|
| 827 | |
|---|
| 828 | !! special MARS |
|---|
| 829 | real r_input(n) |
|---|
| 830 | real cp_input(n) |
|---|
| 831 | real cv_input(n) |
|---|
| 832 | real cvpm_input(n) |
|---|
| 833 | real pfile_input(n) |
|---|
| 834 | real t_input(n) |
|---|
| 835 | real rhofile_input(n) |
|---|
| 836 | !! special MARS |
|---|
| 837 | |
|---|
| 838 | ! diagnostics |
|---|
| 839 | |
|---|
| 840 | real rho_surf, p_input(n), rho_input(n) |
|---|
| 841 | real pm_input(n) ! this are for full moist sounding |
|---|
| 842 | |
|---|
| 843 | ! local data |
|---|
| 844 | |
|---|
| 845 | real p1000mb,cv,cp,r,cvpm,g |
|---|
| 846 | !****Mars |
|---|
| 847 | ! parameter (p1000mb = 1.e+05, r = 287, cp = 1003., cv = cp-r, cvpm = -cv/cp, g=9.81 ) |
|---|
| 848 | ! parameter (p1000mb = 610., r = 192., cp = 844.6, cv = cp-r, cvpm = -cv/cp, g=3.72) |
|---|
| 849 | parameter (p1000mb = 610., r = 191., cp = 744.5, cv = cp-r, cvpm = -cv/cp, g=3.72) |
|---|
| 850 | !****Mars |
|---|
| 851 | integer k, it, nl |
|---|
| 852 | real qvf, qvf1, dz |
|---|
| 853 | |
|---|
| 854 | ! first, read the sounding |
|---|
| 855 | |
|---|
| 856 | call read_sounding( p_surf, th_surf, qv_surf, & |
|---|
| 857 | h_input, th_input, qv_input, u_input, v_input, r_input, cp_input, pfile_input, t_input, rhofile_input, n, nl, debug ) |
|---|
| 858 | |
|---|
| 859 | |
|---|
| 860 | !! special MARS |
|---|
| 861 | do k=1,nl |
|---|
| 862 | cv_input(k) = cp_input(k) - r_input(k) |
|---|
| 863 | cvpm_input(k) = - cv_input(k) / cp_input(k) |
|---|
| 864 | enddo |
|---|
| 865 | !! special MARS |
|---|
| 866 | |
|---|
| 867 | if(dry) then |
|---|
| 868 | do k=1,nl |
|---|
| 869 | qv_input(k) = 0. |
|---|
| 870 | enddo |
|---|
| 871 | endif |
|---|
| 872 | |
|---|
| 873 | if(debug) write(6,*) ' number of input levels = ',nl |
|---|
| 874 | |
|---|
| 875 | nl_in = nl |
|---|
| 876 | if(nl_in .gt. nl_max ) then |
|---|
| 877 | write(6,*) ' too many levels for input arrays ',nl_in,nl_max |
|---|
| 878 | call wrf_error_fatal ( ' too many levels for input arrays ' ) |
|---|
| 879 | end if |
|---|
| 880 | |
|---|
| 881 | ! compute diagnostics, |
|---|
| 882 | ! first, convert qv(g/kg) to qv(g/g) |
|---|
| 883 | |
|---|
| 884 | do k=1,nl |
|---|
| 885 | qv_input(k) = 0.001*qv_input(k) |
|---|
| 886 | enddo |
|---|
| 887 | |
|---|
| 888 | p_surf = 100.*p_surf ! convert to pascals |
|---|
| 889 | qvf = 1. + rvovrd*qv_input(1) |
|---|
| 890 | rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm)) |
|---|
| 891 | pi_surf = (p_surf/p1000mb)**(r/cp) |
|---|
| 892 | !!!!!! rcp variable |
|---|
| 893 | !rho_surf = 1./((r_input(1)/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm_input(1))) |
|---|
| 894 | !pi_surf = (p_surf/p1000mb)**(r_input(1)/cp_input(1)) |
|---|
| 895 | |
|---|
| 896 | |
|---|
| 897 | if(debug) then |
|---|
| 898 | write(6,*) ' surface density is ',rho_surf |
|---|
| 899 | write(6,*) ' surface pi is ',pi_surf |
|---|
| 900 | end if |
|---|
| 901 | |
|---|
| 902 | |
|---|
| 903 | ! integrate moist sounding hydrostatically, starting from the |
|---|
| 904 | ! specified surface pressure |
|---|
| 905 | ! -> first, integrate from surface to lowest level |
|---|
| 906 | |
|---|
| 907 | qvf = 1. + rvovrd*qv_input(1) |
|---|
| 908 | qvf1 = 1. + qv_input(1) |
|---|
| 909 | rho_input(1) = rho_surf |
|---|
| 910 | dz = h_input(1) |
|---|
| 911 | do it=1,10 |
|---|
| 912 | !!MARS MARS |
|---|
| 913 | pm_input(1) = p_surf !& |
|---|
| 914 | ! - dz*(0.25*rho_surf+0.75*rho_input(1))*g*qvf1 !!! BEURK |
|---|
| 915 | ! - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1 !! parce que couche 1 tres proche |
|---|
| 916 | rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm)) |
|---|
| 917 | !!!!!!! rcp variable |
|---|
| 918 | !rho_input(1) = 1./((r_input(1)/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm_input(1))) |
|---|
| 919 | enddo |
|---|
| 920 | |
|---|
| 921 | ! integrate up the column |
|---|
| 922 | |
|---|
| 923 | do k=2,nl |
|---|
| 924 | rho_input(k) = rho_input(k-1) |
|---|
| 925 | dz = h_input(k)-h_input(k-1) |
|---|
| 926 | !!!!!!! rcp variable |
|---|
| 927 | !dz = r_input(k) * t_input(k) * (- p_input(k) + p_input(k-1)) / p_input(k) / g |
|---|
| 928 | !!dz = - cp_input(k) * (- t_input(k) + t_input(k-1)) / g |
|---|
| 929 | qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k))) |
|---|
| 930 | qvf = 1. + rvovrd*qv_input(k) ! qv is in g/kg here |
|---|
| 931 | |
|---|
| 932 | print *, 'input', pfile_input(k), rhofile_input(k) |
|---|
| 933 | |
|---|
| 934 | do it=1,10 !!ou moins??? non. !! trop de rho(k-1) donne une pression trop faible puis crash |
|---|
| 935 | !! mais coeff ci-dessous vont varier la pression calculÃe |
|---|
| 936 | pm_input(k) = pm_input(k-1) & |
|---|
| 937 | - dz*(0.75*rho_input(k)+0.25*rho_input(k-1))*g*qvf1 |
|---|
| 938 | !- 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1 |
|---|
| 939 | rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm)) |
|---|
| 940 | !! |
|---|
| 941 | !! marche pas |
|---|
| 942 | !! |
|---|
| 943 | !pm_input(k) = pm_input(k-1) & |
|---|
| 944 | ! - 0.5*dz*(1./rho_input(k)+1./rho_input(k-1))*g*qvf1 |
|---|
| 945 | !rho_input(k) = (r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm) |
|---|
| 946 | !!!!!!! rcp variable |
|---|
| 947 | !rho_input(k) = 1./((r_input(k)/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm_input(k))) |
|---|
| 948 | !print *, p_input(k), pm_input(k),((r_input(k)/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm_input(k))),k |
|---|
| 949 | print *, it, pm_input(k), rho_input(k), dz |
|---|
| 950 | enddo |
|---|
| 951 | enddo |
|---|
| 952 | |
|---|
| 953 | |
|---|
| 954 | ! we have the moist sounding |
|---|
| 955 | |
|---|
| 956 | ! next, compute the dry sounding using p at the highest level from the |
|---|
| 957 | ! moist sounding and integrating down. |
|---|
| 958 | |
|---|
| 959 | p_input(nl) = pm_input(nl) |
|---|
| 960 | |
|---|
| 961 | do k=nl-1,1,-1 |
|---|
| 962 | dz = h_input(k+1)-h_input(k) |
|---|
| 963 | p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g |
|---|
| 964 | enddo |
|---|
| 965 | |
|---|
| 966 | |
|---|
| 967 | do k=1,nl |
|---|
| 968 | |
|---|
| 969 | zk(k) = h_input(k) |
|---|
| 970 | p(k) = pm_input(k) |
|---|
| 971 | p_dry(k) = p_input(k) |
|---|
| 972 | theta(k) = th_input(k) |
|---|
| 973 | rho(k) = rho_input(k) |
|---|
| 974 | u(k) = mulu*u_input(k) + addu |
|---|
| 975 | v(k) = mulv*v_input(k) + addv |
|---|
| 976 | qv(k) = qv_input(k) |
|---|
| 977 | |
|---|
| 978 | !!!! direct input from file |
|---|
| 979 | write(6,*) '*** DIRECT INPUT FROM FILE ***' |
|---|
| 980 | p(k) = pfile_input(k) |
|---|
| 981 | p_dry(k) = pfile_input(k) |
|---|
| 982 | rho(k) = rhofile_input(k) |
|---|
| 983 | |
|---|
| 984 | enddo |
|---|
| 985 | |
|---|
| 986 | if(debug) then |
|---|
| 987 | write(6,*) ' sounding ' |
|---|
| 988 | write(6,*) ' k height(m) press (Pa) pd(Pa) theta (K) den(kg/m^3) u(m/s) v(m/s) qv(g/g) ' |
|---|
| 989 | do k=1,nl |
|---|
| 990 | 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) |
|---|
| 991 | enddo |
|---|
| 992 | |
|---|
| 993 | end if |
|---|
| 994 | |
|---|
| 995 | end subroutine get_sounding |
|---|
| 996 | |
|---|
| 997 | !------------------------------------------------------- |
|---|
| 998 | |
|---|
| 999 | subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,r,cp,p,t,rho,n,nl,debug ) |
|---|
| 1000 | implicit none |
|---|
| 1001 | integer n,nl |
|---|
| 1002 | real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n),r(n),cp(n),p(n),t(n),rho(n) |
|---|
| 1003 | logical end_of_file |
|---|
| 1004 | logical debug |
|---|
| 1005 | |
|---|
| 1006 | integer k |
|---|
| 1007 | |
|---|
| 1008 | open(unit=10,file='input_sounding',form='formatted',status='old') |
|---|
| 1009 | rewind(10) |
|---|
| 1010 | read(10,*) ps, ts, qvs |
|---|
| 1011 | if(debug) then |
|---|
| 1012 | write(6,*) ' input sounding surface parameters ' |
|---|
| 1013 | write(6,*) ' surface pressure (mb) ',ps |
|---|
| 1014 | write(6,*) ' surface pot. temp (K) ',ts |
|---|
| 1015 | write(6,*) ' surface mixing ratio (g/kg) ',qvs |
|---|
| 1016 | end if |
|---|
| 1017 | |
|---|
| 1018 | end_of_file = .false. |
|---|
| 1019 | k = 0 |
|---|
| 1020 | |
|---|
| 1021 | do while (.not. end_of_file) |
|---|
| 1022 | |
|---|
| 1023 | read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1) |
|---|
| 1024 | k = k+1 |
|---|
| 1025 | if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k) |
|---|
| 1026 | go to 110 |
|---|
| 1027 | 100 end_of_file = .true. |
|---|
| 1028 | 110 continue |
|---|
| 1029 | enddo |
|---|
| 1030 | |
|---|
| 1031 | |
|---|
| 1032 | !!! special MARS |
|---|
| 1033 | open(unit=11,file='input_therm',form='formatted',status='old') |
|---|
| 1034 | rewind(11) |
|---|
| 1035 | end_of_file = .false. |
|---|
| 1036 | k = 0 |
|---|
| 1037 | do while (.not. end_of_file) |
|---|
| 1038 | |
|---|
| 1039 | read(11,*,end=101) r(k+1), cp(k+1), p(k+1), rho(k+1), t(k+1) |
|---|
| 1040 | write(*,*) k, r(k+1), cp(k+1), p(k+1), rho(k+1), t(k+1) |
|---|
| 1041 | k = k+1 |
|---|
| 1042 | go to 112 |
|---|
| 1043 | 101 end_of_file = .true. |
|---|
| 1044 | 112 continue |
|---|
| 1045 | enddo |
|---|
| 1046 | !!! special MARS |
|---|
| 1047 | |
|---|
| 1048 | |
|---|
| 1049 | |
|---|
| 1050 | nl = k |
|---|
| 1051 | |
|---|
| 1052 | close(unit=10,status = 'keep') |
|---|
| 1053 | |
|---|
| 1054 | end subroutine read_sounding |
|---|
| 1055 | |
|---|
| 1056 | END MODULE module_initialize |
|---|
