chgres_cube 1.14.0
Loading...
Searching...
No Matches
atmosphere.F90
Go to the documentation of this file.
1
4
20
21 use esmf
22
31
32 use atm_input_data, only : lev_input, &
45
46 use model_grid, only : target_grid, &
52
57 regional, &
62
69
70 use write_data, only : write_fv3_atm_header_netcdf, &
71 write_fv3_atm_bndy_data_netcdf, &
72 write_fv3_atm_data_netcdf
73
74 use utilities, only : error_handler
75
76 implicit none
77
78 private
79
80 type(esmf_field) :: dzdt_b4adj_target_grid
81 type(esmf_field), allocatable :: tracers_b4adj_target_grid(:)
82 type(esmf_field) :: ps_b4adj_target_grid
83 type(esmf_field) :: pres_target_grid
84 type(esmf_field) :: pres_b4adj_target_grid
85 type(esmf_field) :: temp_b4adj_target_grid
87 type(esmf_field) :: xwind_target_grid
88 type(esmf_field) :: ywind_target_grid
89 type(esmf_field) :: zwind_target_grid
90 type(esmf_field) :: xwind_b4adj_target_grid
91 type(esmf_field) :: ywind_b4adj_target_grid
92 type(esmf_field) :: zwind_b4adj_target_grid
93 type(esmf_field) :: xwind_s_target_grid
94 type(esmf_field) :: ywind_s_target_grid
95 type(esmf_field) :: zwind_s_target_grid
96 type(esmf_field) :: xwind_w_target_grid
97 type(esmf_field) :: ywind_w_target_grid
98 type(esmf_field) :: zwind_w_target_grid
99
100! Fields associated with thompson microphysics climatological tracers.
101
108
109 public :: atmosphere_driver
110 public :: read_vcoord_info
111
112 contains
113
118 subroutine atmosphere_driver(localpet)
119
120 use mpi_f08
121
122 implicit none
123
124 integer, intent(in) :: localpet
125
126 integer :: isrctermprocessing
127 integer :: rc, n
128
129 type(esmf_regridmethod_flag) :: method
130 type(esmf_routehandle) :: regrid_bl
131
132 real(esmf_kind_r8), parameter :: p0=101325.0
133 real(esmf_kind_r8), parameter :: rd = 287.058
134 real(esmf_kind_r8), parameter :: grav = 9.81
135 real(esmf_kind_r8), parameter :: lapse = -6.5e-03
136
137 real(esmf_kind_r8), parameter :: exponent = rd*lapse/grav
138 real(esmf_kind_r8), parameter :: one_over_exponent = 1.0 / exponent
139
140 real(esmf_kind_r8), pointer :: psptr(:,:), tempptr(:,:,:)
141
142!-----------------------------------------------------------------------------------
143! Read atmospheric fields on the input grid.
144!-----------------------------------------------------------------------------------
145
146 call read_input_atm_data(localpet)
147
148!-----------------------------------------------------------------------------------
149! Read vertical coordinate info for target grid.
150!-----------------------------------------------------------------------------------
151
153
154!-----------------------------------------------------------------------------------
155! Create target grid field objects to hold data before vertical adjustment.
156!-----------------------------------------------------------------------------------
157
159
160!-----------------------------------------------------------------------------------
161! Horizontally interpolate. If specified, use weights from file.
162!-----------------------------------------------------------------------------------
163
164 isrctermprocessing = 1
165
166 if (trim(atm_weight_file) /= "NULL") then
167
168 print*,"- CALL FieldSMMStore FOR ATMOSPHERIC FIELDS."
169
170 call esmf_fieldsmmstore(temp_input_grid, &
173 routehandle=regrid_bl, &
174 srctermprocessing=isrctermprocessing, rc=rc)
175 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
176 call error_handler("IN FieldSMMStore", rc)
177
178 else
179
180 print*,"- CALL FieldRegridStore FOR ATMOSPHERIC FIELDS."
181
182 method=esmf_regridmethod_bilinear
183
184 call esmf_fieldregridstore(temp_input_grid, &
186 polemethod=esmf_polemethod_allavg, &
187 srctermprocessing=isrctermprocessing, &
188 routehandle=regrid_bl, &
189 regridmethod=method, rc=rc)
190 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
191 call error_handler("IN FieldRegridStore", rc)
192
193 endif
194
195 print*,"- CALL Field_Regrid FOR TEMPERATURE."
196 call esmf_fieldregrid(temp_input_grid, &
198 routehandle=regrid_bl, &
199 termorderflag=esmf_termorder_srcseq, &
200 rc=rc)
201 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
202 call error_handler("IN FieldRegrid", rc)
203
204 print*,"- CALL Field_Regrid FOR PRESSURE."
205 call esmf_fieldregrid(pres_input_grid, &
207 routehandle=regrid_bl, &
208 termorderflag=esmf_termorder_srcseq, &
209 rc=rc)
210 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
211 call error_handler("IN FieldRegrid", rc)
212
213 do n = 1, num_tracers_input
214 print*,"- CALL Field_Regrid FOR TRACER ", trim(tracers(n))
215 call esmf_fieldregrid(tracers_input_grid(n), &
217 routehandle=regrid_bl, &
218 termorderflag=esmf_termorder_srcseq, rc=rc)
219 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
220 call error_handler("IN FieldRegrid", rc)
221
222 enddo
223
224 print*,"- CALL Field_Regrid FOR VERTICAL VELOCITY."
225 call esmf_fieldregrid(dzdt_input_grid, &
227 routehandle=regrid_bl, &
228 termorderflag=esmf_termorder_srcseq, rc=rc)
229 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
230 call error_handler("IN FieldRegrid", rc)
231
232 nullify(tempptr)
233 print*,"- CALL FieldGet FOR INPUT GRID VERTICAL VEL."
234 call esmf_fieldget(dzdt_input_grid, &
235 farrayptr=tempptr, rc=rc)
236 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
237 call error_handler("IN FieldGet", rc)
238
239 print*, "MIN MAX W INPUT = ", minval(tempptr), maxval(tempptr)
240
241 nullify(tempptr)
242 print*,"- CALL FieldGet FOR VERTICAL VEL B4ADJ."
243 call esmf_fieldget(dzdt_b4adj_target_grid, &
244 farrayptr=tempptr, rc=rc)
245 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
246 call error_handler("IN FieldGet", rc)
247
248 print*, "MIN MAX W B4ADJ = ", minval(tempptr), maxval(tempptr)
249
250 nullify(psptr)
251 print*,"- CALL FieldGet FOR INPUT SURFACE PRESSURE."
252 call esmf_fieldget(ps_input_grid, &
253 farrayptr=psptr, rc=rc)
254 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
255 call error_handler("IN FieldGet", rc)
256
257!------------------------------------------------------------------------------------
258! Assume standard lapse rate when interpolating pressure (per Phil Pegion).
259!------------------------------------------------------------------------------------
260
261 psptr = (psptr/p0)**exponent
262
263 print*,"- CALL Field_Regrid FOR SURFACE PRESSURE."
264 call esmf_fieldregrid(ps_input_grid, &
266 routehandle=regrid_bl, &
267 termorderflag=esmf_termorder_srcseq, &
268 rc=rc)
269 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
270 call error_handler("IN FieldRegrid", rc)
271
272 nullify(psptr)
273 print*,"- CALL FieldGet FOR INPUT SURFACE PRESSURE B4ADJ."
274 call esmf_fieldget(ps_b4adj_target_grid, &
275 farrayptr=psptr, rc=rc)
276 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
277 call error_handler("IN FieldGet", rc)
278
279 psptr = p0 * psptr**one_over_exponent
280
281 print*,"- CALL Field_Regrid FOR TERRAIN."
282 call esmf_fieldregrid(terrain_input_grid, &
284 routehandle=regrid_bl, &
285 termorderflag=esmf_termorder_srcseq, &
286 rc=rc)
287 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
288 call error_handler("IN FieldRegrid", rc)
289
290 print*,"- CALL Field_Regrid FOR x WIND."
291 call esmf_fieldregrid(xwind_input_grid, &
293 routehandle=regrid_bl, &
294 termorderflag=esmf_termorder_srcseq, rc=rc)
295 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
296 call error_handler("IN FieldRegrid", rc)
297
298 print*,"- CALL Field_Regrid FOR y WIND."
299 call esmf_fieldregrid(ywind_input_grid, &
301 routehandle=regrid_bl, &
302 termorderflag=esmf_termorder_srcseq, rc=rc)
303 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
304 call error_handler("IN FieldRegrid", rc)
305
306 print*,"- CALL Field_Regrid FOR z WIND."
307 call esmf_fieldregrid(zwind_input_grid, &
309 routehandle=regrid_bl, &
310 termorderflag=esmf_termorder_srcseq, rc=rc)
311 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
312 call error_handler("IN FieldRegrid", rc)
313
314
315
316 print*,"- CALL FieldRegridRelease."
317 call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
318 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
319 call error_handler("IN FieldRegridRelease", rc)
320
321!-----------------------------------------------------------------------------------
322! Deallocate input fields.
323!-----------------------------------------------------------------------------------
324
326
327!-----------------------------------------------------------------------------------
328! Create target grid field objects to hold data after vertical interpolation.
329!-----------------------------------------------------------------------------------
330
332
333!-----------------------------------------------------------------------------------
334! Adjust surface pressure for terrain differences.
335!-----------------------------------------------------------------------------------
336
337 call newps(localpet)
338
339!-----------------------------------------------------------------------------------
340! Compute 3-d pressure based on adjusted surface pressure.
341!-----------------------------------------------------------------------------------
342
343 call newpr1(localpet)
344
345!-----------------------------------------------------------------------------------
346! Vertically interpolate.
347!-----------------------------------------------------------------------------------
348
349 call vintg
350
351 if( wam_cold_start ) then
353 endif
354
355!-----------------------------------------------------------------------------------
356! Compute height.
357!-----------------------------------------------------------------------------------
358
359 call compute_zh
360
361!-----------------------------------------------------------------------------------
362! Free up memory.
363!-----------------------------------------------------------------------------------
364
366
367!-----------------------------------------------------------------------------------
368! Interpolate winds to 'd' grid.
369!-----------------------------------------------------------------------------------
370
371 isrctermprocessing = 1
372 method=esmf_regridmethod_bilinear
373
374 print*,"- CALL FieldRegridStore FOR X-WIND WEST EDGE."
375 call esmf_fieldregridstore(xwind_target_grid, &
377 polemethod=esmf_polemethod_allavg, &
378 srctermprocessing=isrctermprocessing, &
379 routehandle=regrid_bl, &
380 extrapmethod=esmf_extrapmethod_nearest_stod, &
381 regridmethod=method, rc=rc)
382 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
383 call error_handler("IN FieldRegridStore", rc)
384
385 print*,"- CALL Field_Regrid FOR X-WIND WEST EDGE."
386 call esmf_fieldregrid(xwind_target_grid, &
388 routehandle=regrid_bl, &
389 termorderflag=esmf_termorder_srcseq, rc=rc)
390 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
391 call error_handler("IN FieldRegrid", rc)
392
393 print*,"- CALL Field_Regrid FOR Y-WIND WEST EDGE."
394 call esmf_fieldregrid(ywind_target_grid, &
396 routehandle=regrid_bl, &
397 termorderflag=esmf_termorder_srcseq, rc=rc)
398 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
399 call error_handler("IN FieldRegrid", rc)
400
401 print*,"- CALL Field_Regrid FOR Z-WIND WEST EDGE."
402 call esmf_fieldregrid(zwind_target_grid, &
404 routehandle=regrid_bl, &
405 termorderflag=esmf_termorder_srcseq, rc=rc)
406 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
407 call error_handler("IN FieldRegrid", rc)
408
409 print*,"- CALL FieldRegridRelease."
410 call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
411 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
412 call error_handler("IN FieldRegridRelease", rc)
413
414 isrctermprocessing = 1
415 method=esmf_regridmethod_bilinear
416
417 print*,"- CALL FieldRegridStore FOR X-WIND SOUTH EDGE."
418 call esmf_fieldregridstore(xwind_target_grid, &
420 polemethod=esmf_polemethod_allavg, &
421 srctermprocessing=isrctermprocessing, &
422 routehandle=regrid_bl, &
423 extrapmethod=esmf_extrapmethod_nearest_stod, &
424 regridmethod=method, rc=rc)
425 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
426 call error_handler("IN FieldRegridStore", rc)
427
428 print*,"- CALL Field_Regrid FOR X-WIND SOUTH EDGE."
429 call esmf_fieldregrid(xwind_target_grid, &
431 routehandle=regrid_bl, &
432 termorderflag=esmf_termorder_srcseq, rc=rc)
433 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
434 call error_handler("IN FieldRegrid", rc)
435
436 print*,"- CALL Field_Regrid FOR Y-WIND SOUTH EDGE."
437 call esmf_fieldregrid(ywind_target_grid, &
439 routehandle=regrid_bl, &
440 termorderflag=esmf_termorder_srcseq, rc=rc)
441 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
442 call error_handler("IN FieldRegrid", rc)
443
444 print*,"- CALL Field_Regrid FOR Z-WIND SOUTH EDGE."
445 call esmf_fieldregrid(zwind_target_grid, &
447 routehandle=regrid_bl, &
448 termorderflag=esmf_termorder_srcseq, rc=rc)
449 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
450 call error_handler("IN FieldRegrid", rc)
451
452 print*,"- CALL FieldRegridRelease."
453 call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
454 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
455 call error_handler("IN FieldRegridRelease", rc)
456
457!-----------------------------------------------------------------------------------
458! Convert from 3-d to 2-d cartesian winds.
459!-----------------------------------------------------------------------------------
460
462
463!-----------------------------------------------------------------------------------
464! If selected, process thompson microphysics climatological fields.
465!-----------------------------------------------------------------------------------
466
467 if (use_thomp_mp_climo) then
471 endif
472
473!-----------------------------------------------------------------------------------
474! Write target data to file.
475!-----------------------------------------------------------------------------------
476
477 call write_fv3_atm_header_netcdf(localpet)
478 if (regional <= 1) call write_fv3_atm_data_netcdf(localpet)
479 if (regional >= 1) call write_fv3_atm_bndy_data_netcdf(localpet)
480
481!-----------------------------------------------------------------------------------
482! Free up memory.
483!-----------------------------------------------------------------------------------
484
486
487 end subroutine atmosphere_driver
488
495
496 implicit none
497
498 integer :: rc, n
499
501
502 do n = 1, num_tracers_input
503 print*,"- CALL FieldCreate FOR TARGET GRID TRACER BEFORE ADJUSTMENT ", trim(tracers(n))
504 tracers_b4adj_target_grid(n) = esmf_fieldcreate(target_grid, &
505 typekind=esmf_typekind_r8, &
506 staggerloc=esmf_staggerloc_center, &
507 ungriddedlbound=(/1/), &
508 ungriddedubound=(/lev_input/), rc=rc)
509 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
510 call error_handler("IN FieldCreate", rc)
511 enddo
512
513 print*,"- CALL FieldCreate FOR TARGET GRID TEMPERATURE BEFORE ADJUSTMENT."
514 temp_b4adj_target_grid = esmf_fieldcreate(target_grid, &
515 typekind=esmf_typekind_r8, &
516 staggerloc=esmf_staggerloc_center, &
517 ungriddedlbound=(/1/), &
518 ungriddedubound=(/lev_input/), rc=rc)
519 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
520 call error_handler("IN FieldCreate", rc)
521
522 print*,"- CALL FieldCreate FOR TARGET GRID PRESSURE BEFORE ADJUSTMENT."
523 pres_b4adj_target_grid = esmf_fieldcreate(target_grid, &
524 typekind=esmf_typekind_r8, &
525 staggerloc=esmf_staggerloc_center, &
526 ungriddedlbound=(/1/), &
527 ungriddedubound=(/lev_input/), rc=rc)
528 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
529 call error_handler("IN FieldCreate", rc)
530
531 print*,"- CALL FieldCreate FOR TARGET GRID VERTICAL VELOCITY BEFORE ADJUSTMENT."
532 dzdt_b4adj_target_grid = esmf_fieldcreate(target_grid, &
533 typekind=esmf_typekind_r8, &
534 staggerloc=esmf_staggerloc_center, &
535 ungriddedlbound=(/1/), &
536 ungriddedubound=(/lev_input/), rc=rc)
537 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
538 call error_handler("IN FieldCreate", rc)
539
540 print*,"- CALL FieldCreate FOR TARGET GRID xwind."
541 xwind_b4adj_target_grid = esmf_fieldcreate(target_grid, &
542 typekind=esmf_typekind_r8, &
543 staggerloc=esmf_staggerloc_center, &
544 ungriddedlbound=(/1/), &
545 ungriddedubound=(/lev_input/), rc=rc)
546 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
547 call error_handler("IN FieldCreate", rc)
548
549 print*,"- CALL FieldCreate FOR TARGET GRID ywind."
550 ywind_b4adj_target_grid = esmf_fieldcreate(target_grid, &
551 typekind=esmf_typekind_r8, &
552 staggerloc=esmf_staggerloc_center, &
553 ungriddedlbound=(/1/), &
554 ungriddedubound=(/lev_input/), rc=rc)
555 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
556 call error_handler("IN FieldCreate", rc)
557
558 print*,"- CALL FieldCreate FOR TARGET GRID zwind."
559 zwind_b4adj_target_grid = esmf_fieldcreate(target_grid, &
560 typekind=esmf_typekind_r8, &
561 staggerloc=esmf_staggerloc_center, &
562 ungriddedlbound=(/1/), &
563 ungriddedubound=(/lev_input/), rc=rc)
564 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
565 call error_handler("IN FieldCreate", rc)
566
567 print*,"- CALL FieldCreate FOR TARGET TERRAIN."
568 terrain_interp_to_target_grid = esmf_fieldcreate(target_grid, &
569 typekind=esmf_typekind_r8, &
570 staggerloc=esmf_staggerloc_center, rc=rc)
571 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
572 call error_handler("IN FieldCreate", rc)
573
574 print*,"- CALL FieldCreate FOR TARGET SURFACE PRESSURE BEFORE ADJUSTMENT."
575 ps_b4adj_target_grid = esmf_fieldcreate(target_grid, &
576 typekind=esmf_typekind_r8, &
577 staggerloc=esmf_staggerloc_center, rc=rc)
578 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
579 call error_handler("IN FieldCreate", rc)
580
581 end subroutine create_atm_b4adj_esmf_fields
582
587
588 implicit none
589
590 integer :: rc, n
591
593
594 do n = 1, num_tracers
595 print*,"- CALL FieldCreate FOR TARGET GRID TRACERS ", trim(tracers(n))
596 tracers_target_grid(n) = esmf_fieldcreate(target_grid, &
597 typekind=esmf_typekind_r8, &
598 staggerloc=esmf_staggerloc_center, &
599 ungriddedlbound=(/1/), &
600 ungriddedubound=(/lev_target/), rc=rc)
601 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
602 call error_handler("IN FieldCreate", rc)
603 enddo
604
605 print*,"- CALL FieldCreate FOR TARGET GRID TEMPERATURE."
606 temp_target_grid = esmf_fieldcreate(target_grid, &
607 typekind=esmf_typekind_r8, &
608 staggerloc=esmf_staggerloc_center, &
609 ungriddedlbound=(/1/), &
610 ungriddedubound=(/lev_target/), rc=rc)
611 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
612 call error_handler("IN FieldCreate", rc)
613
614 print*,"- CALL FieldCreate FOR TARGET GRID PRESSURE."
615 pres_target_grid = esmf_fieldcreate(target_grid, &
616 typekind=esmf_typekind_r8, &
617 staggerloc=esmf_staggerloc_center, &
618 ungriddedlbound=(/1/), &
619 ungriddedubound=(/lev_target/), rc=rc)
620 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
621 call error_handler("IN FieldCreate", rc)
622
623 print*,"- CALL FieldCreate FOR TARGET GRID VERTICAL VELOCITY."
624 dzdt_target_grid = esmf_fieldcreate(target_grid, &
625 typekind=esmf_typekind_r8, &
626 staggerloc=esmf_staggerloc_center, &
627 ungriddedlbound=(/1/), &
628 ungriddedubound=(/lev_target/), rc=rc)
629 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
630 call error_handler("IN FieldCreate", rc)
631
632 print*,"- CALL FieldCreate FOR TARGET GRID DELP."
633 delp_target_grid = esmf_fieldcreate(target_grid, &
634 typekind=esmf_typekind_r8, &
635 staggerloc=esmf_staggerloc_center, &
636 ungriddedlbound=(/1/), &
637 ungriddedubound=(/lev_target/), rc=rc)
638 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
639 call error_handler("IN FieldCreate", rc)
640
641 print*,"- CALL FieldCreate FOR TARGET GRID xwind."
642 xwind_target_grid = esmf_fieldcreate(target_grid, &
643 typekind=esmf_typekind_r8, &
644 staggerloc=esmf_staggerloc_center, &
645 ungriddedlbound=(/1/), &
646 ungriddedubound=(/lev_target/), rc=rc)
647 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
648 call error_handler("IN FieldCreate", rc)
649
650 print*,"- CALL FieldCreate FOR TARGET GRID ywind."
651 ywind_target_grid = esmf_fieldcreate(target_grid, &
652 typekind=esmf_typekind_r8, &
653 staggerloc=esmf_staggerloc_center, &
654 ungriddedlbound=(/1/), &
655 ungriddedubound=(/lev_target/), rc=rc)
656 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
657 call error_handler("IN FieldCreate", rc)
658
659 print*,"- CALL FieldCreate FOR TARGET GRID zwind."
660 zwind_target_grid = esmf_fieldcreate(target_grid, &
661 typekind=esmf_typekind_r8, &
662 staggerloc=esmf_staggerloc_center, &
663 ungriddedlbound=(/1/), &
664 ungriddedubound=(/lev_target/), rc=rc)
665 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
666 call error_handler("IN FieldCreate", rc)
667
668 print*,"- CALL FieldCreate FOR TARGET HEIGHT."
669 zh_target_grid = esmf_fieldcreate(target_grid, &
670 typekind=esmf_typekind_r8, &
671 staggerloc=esmf_staggerloc_center, &
672 ungriddedlbound=(/1/), &
673 ungriddedubound=(/levp1_target/), rc=rc)
674 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
675 call error_handler("IN FieldCreate", rc)
676
677 print*,"- CALL FieldCreate FOR TARGET U_S."
678 u_s_target_grid = esmf_fieldcreate(target_grid, &
679 typekind=esmf_typekind_r8, &
680 staggerloc=esmf_staggerloc_edge2, &
681 ungriddedlbound=(/1/), &
682 ungriddedubound=(/lev_target/), rc=rc)
683 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
684 call error_handler("IN FieldCreate", rc)
685
686 print*,"- CALL FieldCreate FOR TARGET V_S."
687 v_s_target_grid = esmf_fieldcreate(target_grid, &
688 typekind=esmf_typekind_r8, &
689 staggerloc=esmf_staggerloc_edge2, &
690 ungriddedlbound=(/1/), &
691 ungriddedubound=(/lev_target/), rc=rc)
692 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
693 call error_handler("IN FieldCreate", rc)
694
695 print*,"- CALL FieldCreate FOR TARGET xwind_S."
696 xwind_s_target_grid = esmf_fieldcreate(target_grid, &
697 typekind=esmf_typekind_r8, &
698 staggerloc=esmf_staggerloc_edge2, &
699 ungriddedlbound=(/1/), &
700 ungriddedubound=(/lev_target/), rc=rc)
701 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
702 call error_handler("IN FieldCreate", rc)
703
704 print*,"- CALL FieldCreate FOR TARGET ywind_S."
705 ywind_s_target_grid = esmf_fieldcreate(target_grid, &
706 typekind=esmf_typekind_r8, &
707 staggerloc=esmf_staggerloc_edge2, &
708 ungriddedlbound=(/1/), &
709 ungriddedubound=(/lev_target/), rc=rc)
710 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
711 call error_handler("IN FieldCreate", rc)
712
713 print*,"- CALL FieldCreate FOR TARGET zwind_S."
714 zwind_s_target_grid = esmf_fieldcreate(target_grid, &
715 typekind=esmf_typekind_r8, &
716 staggerloc=esmf_staggerloc_edge2, &
717 ungriddedlbound=(/1/), &
718 ungriddedubound=(/lev_target/), rc=rc)
719 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
720 call error_handler("IN FieldCreate", rc)
721
722 print*,"- CALL FieldCreate FOR TARGET U_W."
723 u_w_target_grid = esmf_fieldcreate(target_grid, &
724 typekind=esmf_typekind_r8, &
725 staggerloc=esmf_staggerloc_edge1, &
726 ungriddedlbound=(/1/), &
727 ungriddedubound=(/lev_target/), rc=rc)
728 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
729 call error_handler("IN FieldCreate", rc)
730
731 print*,"- CALL FieldCreate FOR TARGET V_W."
732 v_w_target_grid = esmf_fieldcreate(target_grid, &
733 typekind=esmf_typekind_r8, &
734 staggerloc=esmf_staggerloc_edge1, &
735 ungriddedlbound=(/1/), &
736 ungriddedubound=(/lev_target/), rc=rc)
737 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
738 call error_handler("IN FieldCreate", rc)
739
740 print*,"- CALL FieldCreate FOR TARGET xwind_W."
741 xwind_w_target_grid = esmf_fieldcreate(target_grid, &
742 typekind=esmf_typekind_r8, &
743 staggerloc=esmf_staggerloc_edge1, &
744 ungriddedlbound=(/1/), &
745 ungriddedubound=(/lev_target/), rc=rc)
746 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
747 call error_handler("IN FieldCreate", rc)
748
749 print*,"- CALL FieldCreate FOR TARGET ywind_W."
750 ywind_w_target_grid = esmf_fieldcreate(target_grid, &
751 typekind=esmf_typekind_r8, &
752 staggerloc=esmf_staggerloc_edge1, &
753 ungriddedlbound=(/1/), &
754 ungriddedubound=(/lev_target/), rc=rc)
755 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
756 call error_handler("IN FieldCreate", rc)
757
758 print*,"- CALL FieldCreate FOR TARGET zwind_W."
759 zwind_w_target_grid = esmf_fieldcreate(target_grid, &
760 typekind=esmf_typekind_r8, &
761 staggerloc=esmf_staggerloc_edge1, &
762 ungriddedlbound=(/1/), &
763 ungriddedubound=(/lev_target/), rc=rc)
764 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
765 call error_handler("IN FieldCreate", rc)
766
767 print*,"- CALL FieldCreate FOR TARGET SURFACE PRESSURE."
768 ps_target_grid = esmf_fieldcreate(target_grid, &
769 typekind=esmf_typekind_r8, &
770 staggerloc=esmf_staggerloc_center, rc=rc)
771 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
772 call error_handler("IN FieldCreate", rc)
773
774 end subroutine create_atm_esmf_fields
775
780
781 implicit none
782
783 integer :: clb(3), cub(3)
784 integer :: i, j, k, rc
785
786 real(esmf_kind_r8), pointer :: latptr(:,:)
787 real(esmf_kind_r8), pointer :: lonptr(:,:)
788 real(esmf_kind_r8), pointer :: uptr(:,:,:)
789 real(esmf_kind_r8), pointer :: vptr(:,:,:)
790 real(esmf_kind_r8), pointer :: xwindptr(:,:,:)
791 real(esmf_kind_r8), pointer :: ywindptr(:,:,:)
792 real(esmf_kind_r8), pointer :: zwindptr(:,:,:)
793 real(esmf_kind_r8) :: latrad, lonrad
794
795!-----------------------------------------------------------------------------------
796! Convert from 3-d cartesian to 2-cartesian winds
797!-----------------------------------------------------------------------------------
798
799 print*,'- CONVERT WINDS.'
800
801 print*,"- CALL FieldGet FOR xwind_S."
802 call esmf_fieldget(xwind_s_target_grid, &
803 computationallbound=clb, &
804 computationalubound=cub, &
805 farrayptr=xwindptr, rc=rc)
806 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
807 call error_handler("IN FieldGet", rc)
808
809 print*,"- CALL FieldGet FOR ywind_S."
810 call esmf_fieldget(ywind_s_target_grid, &
811 farrayptr=ywindptr, rc=rc)
812 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
813 call error_handler("IN FieldGet", rc)
814
815 print*,"- CALL FieldGet FOR zwind_S."
816 call esmf_fieldget(zwind_s_target_grid, &
817 farrayptr=zwindptr, rc=rc)
818 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
819 call error_handler("IN FieldGet", rc)
820
821 print*,"- CALL FieldGet FOR U_S."
822 call esmf_fieldget(u_s_target_grid, &
823 farrayptr=uptr, rc=rc)
824 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
825 call error_handler("IN FieldGet", rc)
826
827 print*,"- CALL FieldGet FOR V_S."
828 call esmf_fieldget(v_s_target_grid, &
829 farrayptr=vptr, rc=rc)
830 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
831 call error_handler("IN FieldGet", rc)
832
833 print*,"- CALL FieldGet FOR LATITUDE_S."
834 call esmf_fieldget(latitude_s_target_grid, &
835 farrayptr=latptr, rc=rc)
836 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
837 call error_handler("IN FieldGet", rc)
838
839 print*,"- CALL FieldGet FOR LONGITUDE_S."
840 call esmf_fieldget(longitude_s_target_grid, &
841 farrayptr=lonptr, rc=rc)
842 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
843 call error_handler("IN FieldGet", rc)
844
845 do i = clb(1), cub(1)
846 do j = clb(2), cub(2)
847 latrad = latptr(i,j) * acos(-1.) / 180.0
848 lonrad = lonptr(i,j) * acos(-1.) / 180.0
849 do k = clb(3), cub(3)
850 uptr(i,j,k) = xwindptr(i,j,k) * cos(lonrad) + ywindptr(i,j,k) * sin(lonrad)
851 vptr(i,j,k) = -xwindptr(i,j,k) * sin(latrad) * sin(lonrad) + &
852 ywindptr(i,j,k) * sin(latrad) * cos(lonrad) + &
853 zwindptr(i,j,k) * cos(latrad)
854 enddo
855 enddo
856 enddo
857
858 print*,"- CALL FieldGet FOR xwind_w."
859 call esmf_fieldget(xwind_w_target_grid, &
860 computationallbound=clb, &
861 computationalubound=cub, &
862 farrayptr=xwindptr, rc=rc)
863 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
864 call error_handler("IN FieldGet", rc)
865
866 print*,"- CALL FieldGet FOR ywind_w."
867 call esmf_fieldget(ywind_w_target_grid, &
868 farrayptr=ywindptr, rc=rc)
869 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
870 call error_handler("IN FieldGet", rc)
871
872 print*,"- CALL FieldGet FOR zwind_w."
873 call esmf_fieldget(zwind_w_target_grid, &
874 farrayptr=zwindptr, rc=rc)
875 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
876 call error_handler("IN FieldGet", rc)
877
878 print*,"- CALL FieldGet FOR U_W."
879 call esmf_fieldget(u_w_target_grid, &
880 farrayptr=uptr, rc=rc)
881 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
882 call error_handler("IN FieldGet", rc)
883
884 print*,"- CALL FieldGet FOR V_W."
885 call esmf_fieldget(v_w_target_grid, &
886 farrayptr=vptr, rc=rc)
887 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
888 call error_handler("IN FieldGet", rc)
889
890 print*,"- CALL FieldGet FOR LATITUDE_W."
891 call esmf_fieldget(latitude_w_target_grid, &
892 farrayptr=latptr, rc=rc)
893 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
894 call error_handler("IN FieldGet", rc)
895
896 print*,"- CALL FieldGet FOR LONGITUDE_W."
897 call esmf_fieldget(longitude_w_target_grid, &
898 farrayptr=lonptr, rc=rc)
899 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
900 call error_handler("IN FieldGet", rc)
901
902 do i = clb(1), cub(1)
903 do j = clb(2), cub(2)
904 latrad = latptr(i,j) * acos(-1.) / 180.0
905 lonrad = lonptr(i,j) * acos(-1.) / 180.0
906 do k = clb(3), cub(3)
907 uptr(i,j,k) = xwindptr(i,j,k) * cos(lonrad) + ywindptr(i,j,k) * sin(lonrad)
908 vptr(i,j,k) = -xwindptr(i,j,k) * sin(latrad) * sin(lonrad) + &
909 ywindptr(i,j,k) * sin(latrad) * cos(lonrad) + &
910 zwindptr(i,j,k) * cos(latrad)
911 enddo
912 enddo
913 enddo
914
915 end subroutine convert_winds_to_uv
916
948 subroutine newpr1(localpet)
949 implicit none
950
951 integer, intent(in) :: localpet
952
953 integer :: idsl, idvc, rc
954 integer :: i, j, k, clb(3), cub(3)
955
956 real(esmf_kind_r8), parameter :: rd=287.05
957 real(esmf_kind_r8), parameter :: cp=1004.6
958 real(esmf_kind_r8), parameter :: rocp=rd/cp
959 real(esmf_kind_r8), parameter :: rocp1=rocp+1
960 real(esmf_kind_r8), parameter :: rocpr=1/rocp
961
962 real(esmf_kind_r8), pointer :: delp_ptr(:,:,:)
963 real(esmf_kind_r8), pointer :: pptr(:,:,:) ! adjusted 3-d p.
964 real(esmf_kind_r8), pointer :: psptr(:,:) ! adjusted surface p.
965 real(esmf_kind_r8) :: ak, bk
966 real(esmf_kind_r8), allocatable :: pi(:,:,:)
967
968 print*,"COMPUTE 3-D PRESSURE FROM ADJUSTED SURFACE PRESSURE."
969
970 idvc = 2 ! hard wire for now.
971 idsl = 2 ! hard wire for now.
972
973 print*,"- CALL FieldGet FOR 3-D PRES."
974 call esmf_fieldget(pres_target_grid, &
975 computationallbound=clb, &
976 computationalubound=cub, &
977 farrayptr=pptr, rc=rc)
978 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
979 call error_handler("IN FieldGet", rc)
980
981 print*,"- CALL FieldGet FOR DELP."
982 call esmf_fieldget(delp_target_grid, &
983 computationallbound=clb, &
984 computationalubound=cub, &
985 farrayptr=delp_ptr, rc=rc)
986 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
987 call error_handler("IN FieldGet", rc)
988
989 print*,"- CALL FieldGet FOR SURFACE PRESSURE AFTER ADJUSTMENT"
990 call esmf_fieldget(ps_target_grid, &
991 farrayptr=psptr, rc=rc)
992 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
993 call error_handler("IN FieldGet", rc)
994
995 allocate(pi(clb(1):cub(1),clb(2):cub(2),1:levp1_target))
996
997 if(idvc.eq.2) then
998 do k=1,levp1_target
999 ak = vcoord_target(k,1)
1000 bk = vcoord_target(k,2)
1001 do i= clb(1), cub(1)
1002 do j= clb(2), cub(2)
1003 pi(i,j,k) = ak + bk*psptr(i,j)
1004 enddo
1005 enddo
1006 enddo
1007 do k=1,lev_target
1008 do i= clb(1), cub(1)
1009 do j= clb(2), cub(2)
1010 delp_ptr(i,j,k) = pi(i,j,k) - pi(i,j,k+1)
1011 enddo
1012 enddo
1013 enddo
1014 else
1015 call error_handler("PROGRAM ONLY WORKS WITH IDVC 2", 1)
1016 endif
1017
1018 if(idsl.eq.2) then
1019 do k=1,lev_target
1020 do i= clb(1), cub(1)
1021 do j= clb(2), cub(2)
1022 pptr(i,j,k) = (pi(i,j,k)+pi(i,j,k+1))/2.0
1023 enddo
1024 enddo
1025 enddo
1026 else
1027 do k=1,lev_target
1028 do i= clb(1), cub(1)
1029 do j= clb(2), cub(2)
1030 pptr(i,j,k) = ((pi(i,j,k)**rocp1-pi(i,j,k+1)**rocp1)/ &
1031 (rocp1*(pi(i,j,k)-pi(i,j,k+1))))**rocpr
1032 enddo
1033 enddo
1034 enddo
1035 endif
1036
1037 deallocate(pi)
1038
1039 if (localpet == 0) then
1040 print*,'new pres ',pptr(clb(1),clb(2),:)
1041 print*,'delp ',delp_ptr(clb(1),clb(2),:)
1042 endif
1043
1044 end subroutine newpr1
1045
1059 subroutine newps(localpet)
1060
1061 implicit none
1062
1063 integer, intent(in) :: localpet
1064 integer :: i, j, k, ii
1065 integer :: clb(3), cub(3), ls, rc
1066
1067 real(esmf_kind_r8), pointer :: pptr(:,:,:)
1068 real(esmf_kind_r8), pointer :: psptr(:,:)
1069 real(esmf_kind_r8), pointer :: psnewptr(:,:) ! adjusted surface p.
1070 real(esmf_kind_r8), pointer :: tptr(:,:,:)
1071 real(esmf_kind_r8), pointer :: qptr(:,:,:)
1072 real(esmf_kind_r8), pointer :: zsptr(:,:)
1073 real(esmf_kind_r8), pointer :: zsnewptr(:,:)
1074 real(esmf_kind_r8), allocatable :: zu(:,:)
1075 real(esmf_kind_r8), parameter :: beta=-6.5e-3
1076 real(esmf_kind_r8), parameter :: epsilon=1.e-9
1077 real(esmf_kind_r8), parameter :: g=9.80665
1078 real(esmf_kind_r8), parameter :: rd=287.05
1079 real(esmf_kind_r8), parameter :: rv=461.50
1080 real(esmf_kind_r8), parameter :: gor=g/rd
1081 real(esmf_kind_r8), parameter :: fv=rv/rd-1.
1082 real(esmf_kind_r8) :: ftv, fgam, apu, fz0
1083 real(esmf_kind_r8) :: atvu, atv, fz1, fp0
1084 real(esmf_kind_r8) :: apd, azd, agam, azu
1085 real(esmf_kind_r8) :: atvd, fp1, gamma, pu
1086 real(esmf_kind_r8) :: tvu, pd, tvd
1087 real(esmf_kind_r8) :: at, aq, ap, az
1088
1089 ftv(at,aq)=at*(1+fv*aq)
1090 fgam(apu,atvu,apd,atvd)=-gor*log(atvd/atvu)/log(apd/apu)
1091 fz0(ap,atv,azd,apd)=azd+atv/gor*log(apd/ap)
1092 fz1(ap,atv,azd,apd,agam)=azd-atv/agam*((apd/ap)**(-agam/gor)-1)
1093 fp0(az,azu,apu,atvu)=apu*exp(-gor/atvu*(az-azu))
1094 fp1(az,azu,apu,atvu,agam)=apu*(1+agam/atvu*(az-azu))**(-gor/agam)
1095
1096 print*,"- ADJUST SURFACE PRESSURE FOR NEW TERRAIN."
1097
1098 print*,"- CALL FieldGet FOR 3-D PRES."
1099 call esmf_fieldget(pres_b4adj_target_grid, &
1100 computationallbound=clb, &
1101 computationalubound=cub, &
1102 farrayptr=pptr, rc=rc)
1103 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1104 call error_handler("IN FieldGet", rc)
1105
1106 if(localpet==0) then
1107 print*,'old pres ',pptr(clb(1),clb(2),:)
1108 endif
1109
1110 print*,"- CALL FieldGet FOR TEMPERATURE"
1111 call esmf_fieldget(temp_b4adj_target_grid, &
1112 farrayptr=tptr, rc=rc)
1113 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1114 call error_handler("IN FieldGet", rc)
1115
1116! Find specific humidity in the array of tracer fields.
1117
1118 do ii = 1, num_tracers
1119 if (trim(tracers(ii)) == "sphum") exit
1120 enddo
1121
1122 print*,"- CALL FieldGet FOR SPECIFIC HUMIDITY"
1123 call esmf_fieldget(tracers_b4adj_target_grid(ii), &
1124 farrayptr=qptr, rc=rc)
1125 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1126 call error_handler("IN FieldGet", rc)
1127
1128 print*,"- CALL FieldGet FOR SURFACE PRESSURE BEFORE ADJUSTMENT"
1129 call esmf_fieldget(ps_b4adj_target_grid, &
1130 farrayptr=psptr, rc=rc)
1131 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1132 call error_handler("IN FieldGet", rc)
1133
1134 print*,"- CALL FieldGet FOR SURFACE PRESSURE AFTER ADJUSTMENT"
1135 call esmf_fieldget(ps_target_grid, &
1136 farrayptr=psnewptr, rc=rc)
1137 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1138 call error_handler("IN FieldGet", rc)
1139
1140 print*,"- CALL FieldGet FOR OLD TERRAIN"
1141 call esmf_fieldget(terrain_interp_to_target_grid, &
1142 farrayptr=zsptr, rc=rc)
1143 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1144 call error_handler("IN FieldGet", rc)
1145
1146 print*,"- CALL FieldGet FOR NEW TERRAIN"
1147 call esmf_fieldget(terrain_target_grid, &
1148 farrayptr=zsnewptr, rc=rc)
1149 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1150 call error_handler("IN FieldGet", rc)
1151
1152 allocate(zu(clb(1):cub(1),clb(2):cub(2)))
1153
1154!-----------------------------------------------------------------------------------
1155! Note, this routine was adapted from the spectral GFS which labeled the lowest
1156! model layer as '1'.
1157!-----------------------------------------------------------------------------------
1158
1159!-----------------------------------------------------------------------------------
1160! Compute surface pressure below the original ground.
1161!-----------------------------------------------------------------------------------
1162
1163 ls=0
1164 k=1
1165 gamma=beta
1166 do i=clb(1), cub(1)
1167 do j=clb(2), cub(2)
1168 pu=pptr(i,j,k)
1169 tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1170 zu(i,j)=fz1(pu,tvu,zsptr(i,j),psptr(i,j),gamma)
1171 if(zsnewptr(i,j).le.zu(i,j)) then
1172 pu=pptr(i,j,k)
1173 tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1174 if(abs(gamma).gt.epsilon) then
1175 psnewptr(i,j)=fp1(zsnewptr(i,j),zu(i,j),pu,tvu,gamma)
1176 else
1177 psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1178 endif
1179 else
1180 psnewptr(i,j)=0
1181 ls=ls+1
1182 endif
1183 enddo
1184 enddo
1185
1186!-----------------------------------------------------------------------------------
1187! Compute surface pressure above the original ground.
1188!-----------------------------------------------------------------------------------
1189
1190 do k=2,cub(3)
1191 if(ls.gt.0) then
1192 do i=clb(1),cub(1)
1193 do j=clb(2),cub(2)
1194 if(psnewptr(i,j).eq.0) then
1195 pu=pptr(i,j,k)
1196 tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1197 pd=pptr(i,j,k-1)
1198 tvd=ftv(tptr(i,j,k-1),qptr(i,j,k-1))
1199 gamma=fgam(pu,tvu,pd,tvd)
1200 if(abs(gamma).gt.epsilon) then
1201 zu(i,j)=fz1(pu,tvu,zu(i,j),pd,gamma)
1202 else
1203 zu(i,j)=fz0(pu,tvu,zu(i,j),pd)
1204 endif
1205 if(zsnewptr(i,j).le.zu(i,j)) then
1206 if(abs(gamma).gt.epsilon) then
1207 psnewptr(i,j)=fp1(zsnewptr(i,j),zu(i,j),pu,tvu,gamma)
1208 else
1209 psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1210 endif
1211 ls=ls-1
1212 endif
1213 endif
1214 enddo
1215 enddo
1216 endif
1217 enddo
1218
1219!-----------------------------------------------------------------------------------
1220! Compute surface pressure over the top.
1221!-----------------------------------------------------------------------------------
1222
1223
1224 if(ls.gt.0) then
1225 k=cub(3)
1226 gamma=0
1227 do i=clb(1),cub(1)
1228 do j=clb(2),cub(2)
1229 if(psnewptr(i,j).eq.0) then
1230 pu=pptr(i,j,k)
1231 tvu=ftv(tptr(i,j,k),qptr(i,j,k))
1232 psnewptr(i,j)=fp0(zsnewptr(i,j),zu(i,j),pu,tvu)
1233 endif
1234 enddo
1235 enddo
1236 endif
1237
1238 deallocate(zu)
1239
1240 if (localpet == 0) then
1241! do i=clb(1),cub(1)
1242! do j=clb(2),cub(2)
1243 do i=clb(1),clb(1)
1244 do j=clb(2),clb(2)
1245 print*,'sfcp adjust ',(zsnewptr(i,j)-zsptr(i,j)), psptr(i,j),psnewptr(i,j)
1246 enddo
1247 enddo
1248 endif
1249
1250 end subroutine newps
1251
1257 implicit none
1258
1259 integer :: istat, n, k, localpet, idum(2)
1260
1261 real(esmf_kind_r8), allocatable :: dum1d(:)
1262
1263 type(esmf_vm) :: vm
1264
1265 call esmf_vmgetglobal(vm, rc=istat)
1266 if(esmf_logfounderror(rctocheck=istat,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1267 call error_handler("IN VMGetGlobal", istat)
1268
1269 call esmf_vmget(vm, localpet=localpet, rc=istat)
1270 if(esmf_logfounderror(rctocheck=istat,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1271 call error_handler("IN VMGet", istat)
1272
1273 if (localpet == 0) then
1274 print*
1275 print*,"OPEN VERTICAL COORD FILE: ", trim(vcoord_file_target_grid)
1276 open(14, file=trim(vcoord_file_target_grid), form='formatted', iostat=istat, action='read')
1277 if (istat /= 0) then
1278 call error_handler("OPENING VERTICAL COORD FILE", istat)
1279 endif
1280
1281 read(14, *, iostat=istat) idum(1), idum(2)
1282 if (istat /= 0) then
1283 call error_handler("READING VERTICAL COORD FILE", istat)
1284 endif
1285 endif
1286
1287 call esmf_vmbroadcast(vm, idum, 2, 0, rc=istat)
1288 if(esmf_logfounderror(rctocheck=istat,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1289 call error_handler("IN VMBroadcast", istat)
1290
1291 nvcoord_target = idum(1)
1292 lev_target = idum(2)
1293
1295
1296 allocate(dum1d(levp1_target*nvcoord_target)) ! esmf broadcast requires a 1d array.
1298
1299 if (localpet == 0) then
1300 read(14, *, iostat=istat) ((vcoord_target(n,k), k=1,nvcoord_target), n=1,levp1_target)
1301 if (istat /= 0) then
1302 call error_handler("READING VERTICAL COORD FILE", istat)
1303 endif
1304 close(14)
1305 dum1d = reshape(vcoord_target, (/nvcoord_target*levp1_target/))
1306 endif
1307
1308 call esmf_vmbroadcast(vm, dum1d, (nvcoord_target*levp1_target), 0, rc=istat)
1309 if(esmf_logfounderror(rctocheck=istat,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1310 call error_handler("IN VMBroadcast", istat)
1311
1312 vcoord_target = reshape(dum1d, (/levp1_target,nvcoord_target/))
1313
1314 deallocate(dum1d)
1315
1316 end subroutine read_vcoord_info
1317
1323
1324 implicit none
1325
1326 integer :: isrctermprocessing, rc
1327
1328 type(esmf_regridmethod_flag) :: method
1329 type(esmf_routehandle) :: regrid_bl
1330
1331 isrctermprocessing=1
1332
1333 print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNIFA BEFORE ADJUSTMENT."
1334 qnifa_climo_b4adj_target_grid = esmf_fieldcreate(target_grid, &
1335 typekind=esmf_typekind_r8, &
1336 staggerloc=esmf_staggerloc_center, &
1337 ungriddedlbound=(/1/), &
1338 ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1339 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1340 call error_handler("IN FieldCreate", rc)
1341
1342 print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNWFA BEFORE ADJUSTMENT."
1343 qnwfa_climo_b4adj_target_grid = esmf_fieldcreate(target_grid, &
1344 typekind=esmf_typekind_r8, &
1345 staggerloc=esmf_staggerloc_center, &
1346 ungriddedlbound=(/1/), &
1347 ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1348 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1349 call error_handler("IN FieldCreate", rc)
1350
1351 print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO PRESSURE BEFORE ADJUSTMENT."
1353 typekind=esmf_typekind_r8, &
1354 staggerloc=esmf_staggerloc_center, &
1355 ungriddedlbound=(/1/), &
1356 ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
1357 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1358 call error_handler("IN FieldCreate", rc)
1359
1360 print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNIFA."
1361 qnifa_climo_target_grid = esmf_fieldcreate(target_grid, &
1362 typekind=esmf_typekind_r8, &
1363 staggerloc=esmf_staggerloc_center, &
1364 ungriddedlbound=(/1/), &
1365 ungriddedubound=(/lev_target/), rc=rc)
1366 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1367 call error_handler("IN FieldCreate", rc)
1368
1369 print*,"- CALL FieldCreate FOR TARGET GRID THOMP CLIMO QNWFA."
1370 qnwfa_climo_target_grid = esmf_fieldcreate(target_grid, &
1371 typekind=esmf_typekind_r8, &
1372 staggerloc=esmf_staggerloc_center, &
1373 ungriddedlbound=(/1/), &
1374 ungriddedubound=(/lev_target/), rc=rc)
1375 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1376 call error_handler("IN FieldCreate", rc)
1377
1378 print*,"- CALL FieldRegridStore FOR THOMPSON CLIMO FIELDS."
1379
1380 method=esmf_regridmethod_bilinear
1381
1382 call esmf_fieldregridstore(qnifa_climo_input_grid, &
1384 polemethod=esmf_polemethod_allavg, &
1385 srctermprocessing=isrctermprocessing, &
1386 routehandle=regrid_bl, &
1387 regridmethod=method, rc=rc)
1388 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1389 call error_handler("IN FieldRegridStore", rc)
1390
1391 print*,"- CALL Field_Regrid FOR THOMP CLIMO QNIFA."
1392 call esmf_fieldregrid(qnifa_climo_input_grid, &
1394 routehandle=regrid_bl, &
1395 termorderflag=esmf_termorder_srcseq, rc=rc)
1396 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1397 call error_handler("IN FieldRegrid", rc)
1398
1399 print*,"- CALL Field_Regrid FOR THOMP CLIMO QNWFA."
1400 call esmf_fieldregrid(qnwfa_climo_input_grid, &
1402 routehandle=regrid_bl, &
1403 termorderflag=esmf_termorder_srcseq, rc=rc)
1404 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1405 call error_handler("IN FieldRegrid", rc)
1406
1407 print*,"- CALL Field_Regrid FOR THOMP PRESSURE."
1408 call esmf_fieldregrid(thomp_pres_climo_input_grid, &
1410 routehandle=regrid_bl, &
1411 termorderflag=esmf_termorder_srcseq, rc=rc)
1412 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1413 call error_handler("IN FieldRegrid", rc)
1414
1415 print*,"- CALL FieldRegridRelease."
1416 call esmf_fieldregridrelease(routehandle=regrid_bl, rc=rc)
1417 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1418 call error_handler("IN FieldRegridRelease", rc)
1419
1420!-----------------------------------------------------------------------------------
1421! Free up input data memory.
1422!-----------------------------------------------------------------------------------
1423
1425
1426 end subroutine horiz_interp_thomp_mp_climo
1427
1435
1436 implicit none
1437
1438 INTEGER :: CLB(3), CUB(3), RC
1439 INTEGER :: IM, KM1, KM2, NT
1440 INTEGER :: I, J, K
1441
1442 REAL(ESMF_KIND_R8), ALLOCATABLE :: Z1(:,:,:), Z2(:,:,:)
1443 REAL(ESMF_KIND_R8), ALLOCATABLE :: C1(:,:,:,:),C2(:,:,:,:)
1444
1445 REAL(ESMF_KIND_R8), POINTER :: QNIFA1PTR(:,:,:) ! input
1446 REAL(ESMF_KIND_R8), POINTER :: QNIFA2PTR(:,:,:) ! target
1447 REAL(ESMF_KIND_R8), POINTER :: QNWFA1PTR(:,:,:) ! input
1448 REAL(ESMF_KIND_R8), POINTER :: QNWFA2PTR(:,:,:) ! target
1449 REAL(ESMF_KIND_R8), POINTER :: P1PTR(:,:,:) ! input pressure
1450 REAL(ESMF_KIND_R8), POINTER :: P2PTR(:,:,:) ! target pressure
1451
1452 print*,"- VERTICALY INTERPOLATE THOMP MP CLIMO TRACERS."
1453
1454 print*,"- CALL FieldGet FOR 3-D THOMP PRES."
1455 call esmf_fieldget(thomp_pres_climo_b4adj_target_grid, &
1456 computationallbound=clb, &
1457 computationalubound=cub, &
1458 farrayptr=p1ptr, rc=rc)
1459 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1460 call error_handler("IN FieldGet", rc)
1461
1462! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1463! The '1'/'2' arrays hold fields before/after interpolation.
1464! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1465
1466 nt= 2 ! number of thomp tracers
1467
1468 ALLOCATE(z1(clb(1):cub(1),clb(2):cub(2),lev_thomp_mp_climo))
1469 ALLOCATE(z2(clb(1):cub(1),clb(2):cub(2),lev_target))
1470 ALLOCATE(c1(clb(1):cub(1),clb(2):cub(2),lev_thomp_mp_climo,nt))
1471 ALLOCATE(c2(clb(1):cub(1),clb(2):cub(2),lev_target,nt))
1472
1473 z1 = -log(p1ptr)
1474
1475 print*,"- CALL FieldGet FOR 3-D ADJUSTED PRESS"
1476 call esmf_fieldget(pres_target_grid, &
1477 farrayptr=p2ptr, rc=rc)
1478 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1479 call error_handler("IN FieldGet", rc)
1480
1481 z2 = -log(p2ptr)
1482
1483!print*,'pres check 1 ', p1ptr(clb(1),clb(2),:)
1484!print*,'pres check 2 ', p2ptr(clb(1),clb(2),:)
1485
1486 print*,"- CALL FieldGet FOR qnifa before vertical adjustment."
1487 call esmf_fieldget(qnifa_climo_b4adj_target_grid, &
1488 farrayptr=qnifa1ptr, rc=rc)
1489 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1490 call error_handler("IN FieldGet", rc)
1491
1492 c1(:,:,:,1) = qnifa1ptr(:,:,:)
1493
1494 print*,"- CALL FieldGet FOR qnwfa before vertical adjustment."
1495 call esmf_fieldget(qnwfa_climo_b4adj_target_grid, &
1496 farrayptr=qnwfa1ptr, rc=rc)
1497 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1498 call error_handler("IN FieldGet", rc)
1499
1500 c1(:,:,:,2) = qnwfa1ptr(:,:,:)
1501
1502! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1503! PERFORM LAGRANGIAN ONE-DIMENSIONAL INTERPOLATION
1504! THAT IS 4TH-ORDER IN INTERIOR, 2ND-ORDER IN OUTSIDE INTERVALS
1505! AND 1ST-ORDER FOR EXTRAPOLATION.
1506! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1507
1508 im = (cub(1)-clb(1)+1) * (cub(2)-clb(2)+1)
1510 km2= lev_target
1511
1512 CALL terp3(im,1,1,1,1,nt,(im*km1),(im*km2), &
1513 km1,im,im,z1,c1,km2,im,im,z2,c2)
1514
1515 print*,"- CALL FieldGet FOR ADJUSTED climo qnifa."
1516 call esmf_fieldget(qnifa_climo_target_grid, &
1517 farrayptr=qnifa2ptr, rc=rc)
1518 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1519 call error_handler("IN FieldGet", rc)
1520
1521 print*,"- CALL FieldGet FOR ADJUSTED climo qnwfa."
1522 call esmf_fieldget(qnwfa_climo_target_grid, &
1523 farrayptr=qnwfa2ptr, rc=rc)
1524 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1525 call error_handler("IN FieldGet", rc)
1526
1527 DO k=1,lev_target
1528 DO i=clb(1),cub(1)
1529 DO j=clb(2),cub(2)
1530 qnifa2ptr(i,j,k) = c2(i,j,k,1)
1531 qnwfa2ptr(i,j,k) = c2(i,j,k,2)
1532 ENDDO
1533 ENDDO
1534 ENDDO
1535
1536 DEALLOCATE (z1, z2, c1, c2)
1537
1538 call esmf_fielddestroy(qnifa_climo_b4adj_target_grid, rc=rc)
1539 call esmf_fielddestroy(qnwfa_climo_b4adj_target_grid, rc=rc)
1540 call esmf_fielddestroy(thomp_pres_climo_b4adj_target_grid, rc=rc)
1541
1542 END SUBROUTINE vintg_thomp_mp_climo
1543
1544
1558 SUBROUTINE vintg_wam (YEAR,MONTH,DAY,HOUR,PF)
1559
1560 IMPLICIT NONE
1561
1562 include 'mpif.h'
1563
1564 INTEGER, INTENT(IN) :: YEAR,MONTH,DAY,HOUR
1565 CHARACTER(*), INTENT(IN) :: PF
1566
1567 REAL(ESMF_KIND_R8), PARAMETER :: AMO = 15.9994 ! molecular weight of o
1568 REAL(ESMF_KIND_R8), PARAMETER :: AMO2 = 31.999 !molecular weight of o2
1569 REAL(ESMF_KIND_R8), PARAMETER :: AMN2 = 28.013 !molecular weight of n2
1570
1571 REAL(ESMF_KIND_R8) :: COE,WFUN(10),DEGLAT,HOLD
1572 REAL(ESMF_KIND_R8) :: SUMMASS,QVMASS,O3MASS
1573 INTEGER :: I, J, K, II, CLB(3), CUB(3), RC, KREF
1574 INTEGER :: IDAT(8),JDOW,JDAY,ICDAY
1575
1576 REAL(ESMF_KIND_R8), ALLOCATABLE :: TEMP(:),ON(:),O2N(:),N2N(:),PRMB(:)
1577
1578 REAL(ESMF_KIND_R8), POINTER :: LATPTR(:,:) ! output latitude
1579 REAL(ESMF_KIND_R8), POINTER :: P1PTR(:,:,:) ! input pressure
1580 REAL(ESMF_KIND_R8), POINTER :: P2PTR(:,:,:) ! output pressure
1581 REAL(ESMF_KIND_R8), POINTER :: DZDT2PTR(:,:,:) ! output vvel
1582 REAL(ESMF_KIND_R8), POINTER :: T2PTR(:,:,:) ! output temperature
1583 REAL(ESMF_KIND_R8), POINTER :: Q2PTR(:,:,:) ! output tracer
1584 REAL(ESMF_KIND_R8), POINTER :: QVPTR(:,:,:) ! output tracer
1585 REAL(ESMF_KIND_R8), POINTER :: QOPTR(:,:,:) ! output tracer
1586 REAL(ESMF_KIND_R8), POINTER :: O2PTR(:,:,:) ! output tracer
1587 REAL(ESMF_KIND_R8), POINTER :: O3PTR(:,:,:) ! output tracer
1588 REAL(ESMF_KIND_R8), POINTER :: XWIND2PTR(:,:,:) ! output wind (x component)
1589 REAL(ESMF_KIND_R8), POINTER :: YWIND2PTR(:,:,:) ! output wind (y component)
1590 REAL(ESMF_KIND_R8), POINTER :: ZWIND2PTR(:,:,:) ! output wind (z component)
1591
1592! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1593
1594 print*,"VINTG_WAM:- VERTICALY EXTEND FIELDS FOR WAM COLD START."
1595
1596! prepare date
1597 idat = 0
1598 jdow = 0
1599 jday = 0
1600 icday = 0
1601 idat(1)=year
1602 idat(2)=month
1603 idat(3)=day
1604 idat(5)=hour
1605 CALL w3doxdat(idat,jdow,icday,jday)
1606 print *,"VINTG_WAM: WAM START DATE FOR ICDAY=",icday
1607
1608! prepare weighting function
1609 DO k=1,10
1610 wfun(k) = (k-1.0) / 9.0
1611 ENDDO
1612
1613 ALLOCATE(temp(lev_target))
1614 ALLOCATE(prmb(lev_target))
1615 ALLOCATE( on(lev_target))
1616 ALLOCATE( o2n(lev_target))
1617 ALLOCATE( n2n(lev_target))
1618
1619! p1 (pascal)
1620 print*,"VINTG_WAM:- CALL FieldGet FOR 3-D PRES."
1621 call esmf_fieldget(pres_b4adj_target_grid, &
1622 computationallbound=clb, &
1623 computationalubound=cub, &
1624 farrayptr=p1ptr, rc=rc)
1625 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1626 call error_handler("IN FieldGet", rc)
1627!print*,"VINTG_WAM: p1ptr ",(p1ptr(1,1,k),k=1,LEV_INPUT)
1628
1629! p2 (pascal)
1630 print*,"VINTG_WAM:- CALL FieldGet FOR 3-D ADJUSTED PRESS"
1631 call esmf_fieldget(pres_target_grid, &
1632 farrayptr=p2ptr, rc=rc)
1633 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1634 call error_handler("IN FieldGet", rc)
1635!print*,"VINTG_WAM: p2ptr ",(p2ptr(1,1,k),k=1,LEV_TARGET)
1636
1637! latitude in degree
1638 print*,"VINTG_WAM - CALL FieldGet FOR LATITUDE_S."
1639 call esmf_fieldget(latitude_s_target_grid, &
1640 farrayptr=latptr, rc=rc)
1641 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1642 call error_handler("IN FieldGet", rc)
1643!print*,"VINTG_WAM: latptr ",(latptr(1,j),j=clb(2),cub(2))
1644
1645! temp
1646 print*,"VINTG_WAM:- CALL FieldGet FOR 3-D ADJUSTED TEMP."
1647 call esmf_fieldget(temp_target_grid, &
1648 farrayptr=t2ptr, rc=rc)
1649 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1650 call error_handler("IN FieldGet", rc)
1651
1652! dzdt
1653 print*,"VINTG_WAM:- CALL FieldGet FOR ADJUSTED VERTICAL VELOCITY."
1654 call esmf_fieldget(dzdt_target_grid, &
1655 farrayptr=dzdt2ptr, rc=rc)
1656 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1657 call error_handler("IN FieldGet", rc)
1658
1659! wind
1660 print*,"VINTG_WAM:- CALL FieldGet FOR ADJUSTED WIND COMPONENTS."
1661
1662 call esmf_fieldget(xwind_target_grid, &
1663 farrayptr=xwind2ptr, rc=rc)
1664 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1665 call error_handler("IN FieldGet", rc)
1666
1667 call esmf_fieldget(ywind_target_grid, &
1668 farrayptr=ywind2ptr, rc=rc)
1669 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1670 call error_handler("IN FieldGet", rc)
1671
1672 call esmf_fieldget(zwind_target_grid, &
1673 farrayptr=zwind2ptr, rc=rc)
1674 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1675 call error_handler("IN FieldGet", rc)
1676
1677!
1678! determine vertical blending point and modified extrapolation values
1679!
1680 DO i=clb(1),cub(1)
1681 DO j=clb(2),cub(2)
1682
1683 DO k=1,lev_target
1684 IF(p2ptr(i,j,k).le.p1ptr(i,j,lev_input)) THEN
1685 kref =k-1
1686 EXIT
1687 ENDIF
1688 ENDDO
1689!
1690 DO k=kref,lev_target
1691 coe = p2ptr(i,j,k) / p2ptr(i,j,kref)
1692 xwind2ptr(i,j,k) = coe*xwind2ptr(i,j,k)
1693 ywind2ptr(i,j,k) = coe*ywind2ptr(i,j,k)
1694 zwind2ptr(i,j,k) = coe*zwind2ptr(i,j,k)
1695 dzdt2ptr(i,j,k) = coe*dzdt2ptr(i,j,k)
1696 ENDDO
1697
1698 ENDDO
1699 ENDDO
1700
1701!
1702! point necessary tracers
1703!
1704 DO ii = 1, num_tracers
1705
1706 print*,"VINTG_WAM:- CALL FieldGet FOR 3-D TRACER ", trim(tracers(ii))
1707 call esmf_fieldget(tracers_target_grid(ii), &
1708 farrayptr=q2ptr, rc=rc)
1709 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1710 call error_handler("IN FieldGet", rc)
1711
1712 DO j=clb(2),cub(2)
1713 DO i=clb(1),cub(1)
1714 DO k=1,lev_target
1715 IF(p2ptr(i,j,k).le.p1ptr(i,j,lev_input)) THEN
1716 kref =k-1
1717 EXIT
1718 ENDIF
1719 ENDDO
1720!
1721 DO k=kref,lev_target
1722 coe = min(1.0, p2ptr(i,j,k) / p2ptr(i,j,kref) )
1723 q2ptr(i,j,k) = coe * q2ptr(i,j,k)
1724 ENDDO
1725 ENDDO
1726 ENDDO
1727
1728 IF (trim(tracers(ii)) == "sphum") qvptr => q2ptr
1729 IF (trim(tracers(ii)) == "spo" ) qoptr => q2ptr
1730 IF (trim(tracers(ii)) == "spo2" ) o2ptr => q2ptr
1731 IF (trim(tracers(ii)) == "spo3" ) o3ptr => q2ptr
1732
1733 ENDDO
1734
1735!
1736! obtained wam gases distribution and temperature profile
1737!
1738 DO i=clb(1),cub(1)
1739 DO j=clb(2),cub(2)
1740!
1741 deglat = latptr(i,j)
1742 DO k=1,lev_target
1743 prmb(k) = p2ptr(i,j,k) * 0.01
1744 ENDDO
1745 CALL gettemp(icday,deglat,prmb,lev_target,pf,temp,on,o2n,n2n)
1746!
1747 DO k=1,lev_target
1748 summass = on(k)*amo+o2n(k)*amo2+n2n(k)*amn2
1749 qvmass = summass*qvptr(i,j,k)/(1.-qvptr(i,j,k))
1750 summass = summass+qvmass
1751 o3mass = summass*o3ptr(i,j,k)
1752 summass = summass+o3mass
1753 hold = 1.0 / summass
1754 qoptr(i,j,k) = on(k)*amo *hold
1755 o2ptr(i,j,k) = o2n(k)*amo2*hold
1756 o3ptr(i,j,k) = o3mass * hold
1757 qvptr(i,j,k) = qvmass * hold
1758 ENDDO
1759!
1760 DO k=1,lev_target
1761 IF(p2ptr(i,j,k).le.p1ptr(i,j,lev_input)) THEN
1762 kref =k-1
1763 EXIT
1764 ENDIF
1765 ENDDO
1766!
1767 DO k=kref,lev_target
1768 t2ptr(i,j,k) = temp(k)
1769 ENDDO
1770 DO k=kref-10,kref-1
1771 t2ptr(i,j,k) = wfun(k-kref+11) * temp(k) + &
1772 (1.- wfun(k-kref+11)) * t2ptr(i,j,k)
1773 ENDDO
1774 ENDDO
1775 ENDDO
1776
1777 DEALLOCATE (temp, prmb, on, o2n, n2n)
1778
1779 END SUBROUTINE vintg_wam
1780
1794 SUBROUTINE vintg
1795 use mpi_f08
1796
1797 IMPLICIT NONE
1798
1799 REAL(ESMF_KIND_R8), PARAMETER :: DLTDZ=-6.5e-3*287.05/9.80665
1800 REAL(ESMF_KIND_R8), PARAMETER :: DLPVDRT=-2.5e6/461.50
1801 REAL(ESMF_KIND_R8), PARAMETER :: ONE = 1.0_esmf_kind_r8
1802
1803 INTEGER :: I, J, K, CLB(3), CUB(3), RC
1804 INTEGER :: IM, KM1, KM2, NT, II
1805
1806 REAL(ESMF_KIND_R8) :: DZ
1807 REAL(ESMF_KIND_R8), ALLOCATABLE :: Z1(:,:,:), Z2(:,:,:)
1808 REAL(ESMF_KIND_R8), ALLOCATABLE :: C1(:,:,:,:),C2(:,:,:,:)
1809
1810 REAL(ESMF_KIND_R8), POINTER :: P1PTR(:,:,:) ! input pressure
1811 REAL(ESMF_KIND_R8), POINTER :: P2PTR(:,:,:) ! output pressure
1812 REAL(ESMF_KIND_R8), POINTER :: DZDT1PTR(:,:,:) ! input vvel
1813 REAL(ESMF_KIND_R8), POINTER :: DZDT2PTR(:,:,:) ! output vvel
1814 REAL(ESMF_KIND_R8), POINTER :: T1PTR(:,:,:) ! input temperature
1815 REAL(ESMF_KIND_R8), POINTER :: T2PTR(:,:,:) ! output temperature
1816 REAL(ESMF_KIND_R8), POINTER :: Q1PTR(:,:,:) ! input tracer
1817 REAL(ESMF_KIND_R8), POINTER :: Q2PTR(:,:,:) ! output tracer
1818 REAL(ESMF_KIND_R8), POINTER :: XWIND1PTR(:,:,:) ! input wind (x component)
1819 REAL(ESMF_KIND_R8), POINTER :: YWIND1PTR(:,:,:) ! input wind (y component)
1820 REAL(ESMF_KIND_R8), POINTER :: ZWIND1PTR(:,:,:) ! input wind (z component)
1821 REAL(ESMF_KIND_R8), POINTER :: XWIND2PTR(:,:,:) ! output wind (x component)
1822 REAL(ESMF_KIND_R8), POINTER :: YWIND2PTR(:,:,:) ! output wind (y component)
1823 REAL(ESMF_KIND_R8), POINTER :: ZWIND2PTR(:,:,:) ! output wind (z component)
1824
1825! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1826! COMPUTE LOG PRESSURE INTERPOLATING COORDINATE
1827! AND COPY INPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS
1828! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1829
1830 print*,"- VERTICALY INTERPOLATE FIELDS."
1831
1832 print*,"- CALL FieldGet FOR 3-D PRES."
1833 call esmf_fieldget(pres_b4adj_target_grid, &
1834 computationallbound=clb, &
1835 computationalubound=cub, &
1836 farrayptr=p1ptr, rc=rc)
1837 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1838 call error_handler("IN FieldGet", rc)
1839
1840! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1841! The '1'/'2' arrays hold fields before/after interpolation.
1842! Note the 'z' component of the horizontal wind will be treated as a
1843! tracer. So add one extra third dimension to these 3-d arrays.
1844! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1845
1846 ALLOCATE(z1(clb(1):cub(1),clb(2):cub(2),lev_input))
1847 ALLOCATE(z2(clb(1):cub(1),clb(2):cub(2),lev_target))
1848 ALLOCATE(c1(clb(1):cub(1),clb(2):cub(2),lev_input,num_tracers_input+5))
1849 ALLOCATE(c2(clb(1):cub(1),clb(2):cub(2),lev_target,num_tracers_input+5))
1850
1851 z1 = -log(p1ptr)
1852
1853 print*,"- CALL FieldGet FOR 3-D ADJUSTED PRESS"
1854 call esmf_fieldget(pres_target_grid, &
1855 farrayptr=p2ptr, rc=rc)
1856 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1857 call error_handler("IN FieldGet", rc)
1858
1859 z2 = -log(p2ptr)
1860
1861 print*,"- CALL FieldGet FOR x WIND."
1862 call esmf_fieldget(xwind_b4adj_target_grid, &
1863 farrayptr=xwind1ptr, rc=rc)
1864 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1865 call error_handler("IN FieldGet", rc)
1866
1867 c1(:,:,:,1) = xwind1ptr(:,:,:)
1868
1869 print*,"- CALL FieldGet FOR y WIND."
1870 call esmf_fieldget(ywind_b4adj_target_grid, &
1871 farrayptr=ywind1ptr, rc=rc)
1872 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1873 call error_handler("IN FieldGet", rc)
1874
1875 c1(:,:,:,2) = ywind1ptr(:,:,:)
1876
1877 print*,"- CALL FieldGet FOR z WIND."
1878 call esmf_fieldget(zwind_b4adj_target_grid, &
1879 farrayptr=zwind1ptr, rc=rc)
1880 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1881 call error_handler("IN FieldGet", rc)
1882
1883 c1(:,:,:,3) = zwind1ptr(:,:,:)
1884
1885 print*,"- CALL FieldGet FOR VERTICAL VELOCITY."
1886 call esmf_fieldget(dzdt_b4adj_target_grid, &
1887 farrayptr=dzdt1ptr, rc=rc)
1888 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1889 call error_handler("IN FieldGet", rc)
1890
1891 c1(:,:,:,4) = dzdt1ptr(:,:,:)
1892 print*,"MIN MAX W TARGETB4 IN VINTG = ", minval(dzdt1ptr(:,:,:)), maxval(dzdt1ptr(:,:,:))
1893
1894 print*,"- CALL FieldGet FOR 3-D TEMP."
1895 call esmf_fieldget(temp_b4adj_target_grid, &
1896 farrayptr=t1ptr, rc=rc)
1897 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1898 call error_handler("IN FieldGet", rc)
1899
1900 c1(:,:,:,5) = t1ptr(:,:,:)
1901
1902 DO i = 1, num_tracers_input
1903
1904 print*,"- CALL FieldGet FOR 3-D TRACERS ", trim(tracers(i))
1905 call esmf_fieldget(tracers_b4adj_target_grid(i), &
1906 farrayptr=q1ptr, rc=rc)
1907 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1908 call error_handler("IN FieldGet", rc)
1909
1910 c1(:,:,:,5+i) = q1ptr(:,:,:)
1911
1912 ENDDO
1913
1914! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1915! PERFORM LAGRANGIAN ONE-DIMENSIONAL INTERPOLATION
1916! THAT IS 4TH-ORDER IN INTERIOR, 2ND-ORDER IN OUTSIDE INTERVALS
1917! AND 1ST-ORDER FOR EXTRAPOLATION.
1918! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1919
1920 im = (cub(1)-clb(1)+1) * (cub(2)-clb(2)+1)
1921 km1= lev_input
1922 km2= lev_target
1923 nt= num_tracers_input + 1 ! treat 'z' wind as tracer.
1924
1925 CALL terp3(im,1,1,1,1,4+nt,(im*km1),(im*km2), &
1926 km1,im,im,z1,c1,km2,im,im,z2,c2)
1927! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1928! COPY OUTPUT WIND, TEMPERATURE, HUMIDITY AND OTHER TRACERS
1929! EXCEPT BELOW THE INPUT DOMAIN, LET TEMPERATURE INCREASE WITH A FIXED
1930! LAPSE RATE AND LET THE RELATIVE HUMIDITY REMAIN CONSTANT.
1931! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1932
1933 print*,"- CALL FieldGet FOR 3-D ADJUSTED TEMP."
1934 call esmf_fieldget(temp_target_grid, &
1935 farrayptr=t2ptr, rc=rc)
1936 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1937 call error_handler("IN FieldGet", rc)
1938
1939 print*,"- CALL FieldGet FOR ADJUSTED VERTICAL VELOCITY."
1940 call esmf_fieldget(dzdt_target_grid, &
1941 farrayptr=dzdt2ptr, rc=rc)
1942 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1943 call error_handler("IN FieldGet", rc)
1944
1945 print*,"- CALL FieldGet FOR ADJUSTED xwind."
1946 call esmf_fieldget(xwind_target_grid, &
1947 farrayptr=xwind2ptr, rc=rc)
1948 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1949 call error_handler("IN FieldGet", rc)
1950
1951 print*,"- CALL FieldGet FOR ADJUSTED ywind."
1952 call esmf_fieldget(ywind_target_grid, &
1953 farrayptr=ywind2ptr, rc=rc)
1954 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1955 call error_handler("IN FieldGet", rc)
1956
1957 print*,"- CALL FieldGet FOR ADJUSTED zwind."
1958 call esmf_fieldget(zwind_target_grid, &
1959 farrayptr=zwind2ptr, rc=rc)
1960 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1961 call error_handler("IN FieldGet", rc)
1962
1963 DO k=1,lev_target
1964 DO i=clb(1),cub(1)
1965 DO j=clb(2),cub(2)
1966 xwind2ptr(i,j,k)=c2(i,j,k,1)
1967 ywind2ptr(i,j,k)=c2(i,j,k,2)
1968 zwind2ptr(i,j,k)=c2(i,j,k,3)
1969 dzdt2ptr(i,j,k)=c2(i,j,k,4)
1970 dz=z2(i,j,k)-z1(i,j,1)
1971 IF(dz.GE.0) THEN
1972 t2ptr(i,j,k)=c2(i,j,k,5)
1973 ELSE
1974 t2ptr(i,j,k)=c1(i,j,1,5)*exp(dltdz*dz)
1975 ENDIF
1976 ENDDO
1977 ENDDO
1978 ENDDO
1979
1980 DO ii = 1, num_tracers_input
1981
1982 print*,"- CALL FieldGet FOR 3-D TRACER ", trim(tracers(ii))
1983 call esmf_fieldget(tracers_target_grid(ii), &
1984 farrayptr=q2ptr, rc=rc)
1985 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1986 call error_handler("IN FieldGet", rc)
1987
1988 IF (trim(tracers(ii)) == "sphum") THEN ! specific humidity
1989
1990 DO k=1,lev_target
1991 DO i=clb(1),cub(1)
1992 DO j=clb(2),cub(2)
1993 dz=z2(i,j,k)-z1(i,j,1)
1994 IF(dz.GE.0) THEN
1995 q2ptr(i,j,k) = c2(i,j,k,5+ii)
1996 ELSE
1997 q2ptr(i,j,k) = c1(i,j,1,5+ii)*exp(dlpvdrt*(one/t2ptr(i,j,k)-one/t1ptr(i,j,1))-dz)
1998 ENDIF
1999 ENDDO
2000 ENDDO
2001 ENDDO
2002
2003 ELSE ! all other tracers
2004
2005 DO k=1,lev_target
2006 DO i=clb(1),cub(1)
2007 DO j=clb(2),cub(2)
2008 q2ptr(i,j,k) = c2(i,j,k,5+ii)
2009 ENDDO
2010 ENDDO
2011 ENDDO
2012
2013 ENDIF
2014
2015 ENDDO
2016
2017 DEALLOCATE (z1, z2, c1, c2)
2018
2019 END SUBROUTINE vintg
2020
2057 SUBROUTINE terp3(IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2, &
2058 KM1,KXZ1,KXQ1,Z1,Q1,KM2,KXZ2,KXQ2,Z2,Q2)
2059 IMPLICIT NONE
2060 INTEGER IM,IXZ1,IXQ1,IXZ2,IXQ2,NM,NXQ1,NXQ2
2061 INTEGER KM1,KXZ1,KXQ1,KM2,KXZ2,KXQ2
2062 INTEGER I,K1,K2,N
2063 INTEGER K1S(IM,KM2)
2064 REAL(ESMF_KIND_R8), PARAMETER :: ONE = 1.0_esmf_kind_r8
2065 REAL(ESMF_KIND_R8) :: Z1(1+(IM-1)*IXZ1+(KM1-1)*KXZ1)
2066 REAL(ESMF_KIND_R8) :: Q1(1+(IM-1)*IXQ1+(KM1-1)*KXQ1+(NM-1)*NXQ1)
2067 REAL(ESMF_KIND_R8) :: Z2(1+(IM-1)*IXZ2+(KM2-1)*KXZ2)
2068 REAL(ESMF_KIND_R8) :: Q2(1+(IM-1)*IXQ2+(KM2-1)*KXQ2+(NM-1)*NXQ2)
2069! REAL(ESMF_KIND_R8) :: J2(1+(IM-1)*IXQ2+(KM2-1)*KXQ2+(NM-1)*NXQ2)
2070 REAL(ESMF_KIND_R8) :: FFA(IM),FFB(IM),FFC(IM),FFD(IM)
2071 REAL(ESMF_KIND_R8) :: GGA(IM),GGB(IM),GGC(IM),GGD(IM)
2072 REAL(ESMF_KIND_R8) :: Z1A,Z1B,Z1C,Z1D,Q1A,Q1B,Q1C,Q1D,Z2S,Q2S
2073! REAL(ESMF_KIND_R8) :: J2S
2074
2075! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2076! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT.
2077 CALL rsearch(im,km1,ixz1,kxz1,z1,km2,ixz2,kxz2,z2,1,im,k1s)
2078
2079! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2080! GENERALLY INTERPOLATE CUBICALLY WITH MONOTONIC CONSTRAINT
2081! FROM TWO NEAREST INPUT POINTS ON EITHER SIDE OF THE OUTPUT POINT,
2082! BUT WITHIN THE TWO EDGE INTERVALS INTERPOLATE LINEARLY.
2083! KEEP THE OUTPUT FIELDS CONSTANT OUTSIDE THE INPUT DOMAIN.
2084
2085!$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(IM,IXZ1,IXQ1,IXZ2), &
2086!$OMP& SHARED(IXQ2,NM,NXQ1,NXQ2,KM1,KXZ1,KXQ1,Z1,Q1,KM2,KXZ2), &
2087!$OMP& SHARED(KXQ2,Z2,Q2,K1S)
2088 DO k2=1,km2
2089 DO i=1,im
2090 k1=k1s(i,k2)
2091 IF(k1.EQ.1.OR.k1.EQ.km1-1) THEN
2092 z2s=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
2093 z1a=z1(1+(i-1)*ixz1+(k1-1)*kxz1)
2094 z1b=z1(1+(i-1)*ixz1+(k1+0)*kxz1)
2095 ffa(i)=(z2s-z1b)/(z1a-z1b)
2096 ffb(i)=(z2s-z1a)/(z1b-z1a)
2097 gga(i)=one/(z1a-z1b)
2098 ggb(i)=one/(z1b-z1a)
2099 ELSEIF(k1.GT.1.AND.k1.LT.km1-1) THEN
2100 z2s=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
2101 z1a=z1(1+(i-1)*ixz1+(k1-2)*kxz1)
2102 z1b=z1(1+(i-1)*ixz1+(k1-1)*kxz1)
2103 z1c=z1(1+(i-1)*ixz1+(k1+0)*kxz1)
2104 z1d=z1(1+(i-1)*ixz1+(k1+1)*kxz1)
2105 ffa(i)=(z2s-z1b)/(z1a-z1b)* &
2106 (z2s-z1c)/(z1a-z1c)* &
2107 (z2s-z1d)/(z1a-z1d)
2108 ffb(i)=(z2s-z1a)/(z1b-z1a)* &
2109 (z2s-z1c)/(z1b-z1c)* &
2110 (z2s-z1d)/(z1b-z1d)
2111 ffc(i)=(z2s-z1a)/(z1c-z1a)* &
2112 (z2s-z1b)/(z1c-z1b)* &
2113 (z2s-z1d)/(z1c-z1d)
2114 ffd(i)=(z2s-z1a)/(z1d-z1a)* &
2115 (z2s-z1b)/(z1d-z1b)* &
2116 (z2s-z1c)/(z1d-z1c)
2117 gga(i)= one/(z1a-z1b)* &
2118 (z2s-z1c)/(z1a-z1c)* &
2119 (z2s-z1d)/(z1a-z1d)+ &
2120 (z2s-z1b)/(z1a-z1b)* &
2121 one/(z1a-z1c)* &
2122 (z2s-z1d)/(z1a-z1d)+ &
2123 (z2s-z1b)/(z1a-z1b)* &
2124 (z2s-z1c)/(z1a-z1c)* &
2125 one/(z1a-z1d)
2126 ggb(i)= one/(z1b-z1a)* &
2127 (z2s-z1c)/(z1b-z1c)* &
2128 (z2s-z1d)/(z1b-z1d)+ &
2129 (z2s-z1a)/(z1b-z1a)* &
2130 one/(z1b-z1c)* &
2131 (z2s-z1d)/(z1b-z1d)+ &
2132 (z2s-z1a)/(z1b-z1a)* &
2133 (z2s-z1c)/(z1b-z1c)* &
2134 one/(z1b-z1d)
2135 ggc(i)= one/(z1c-z1a)* &
2136 (z2s-z1b)/(z1c-z1b)* &
2137 (z2s-z1d)/(z1c-z1d)+ &
2138 (z2s-z1a)/(z1c-z1a)* &
2139 one/(z1c-z1b)* &
2140 (z2s-z1d)/(z1c-z1d)+ &
2141 (z2s-z1a)/(z1c-z1a)* &
2142 (z2s-z1b)/(z1c-z1b)* &
2143 one/(z1c-z1d)
2144 ggd(i)= one/(z1d-z1a)* &
2145 (z2s-z1b)/(z1d-z1b)* &
2146 (z2s-z1c)/(z1d-z1c)+ &
2147 (z2s-z1a)/(z1d-z1a)* &
2148 one/(z1d-z1b)* &
2149 (z2s-z1c)/(z1d-z1c)+ &
2150 (z2s-z1a)/(z1d-z1a)* &
2151 (z2s-z1b)/(z1d-z1b)* &
2152 one/(z1d-z1c)
2153 ENDIF
2154 ENDDO
2155
2156! INTERPOLATE.
2157 DO n=1,nm
2158 DO i=1,im
2159 k1=k1s(i,k2)
2160 IF(k1.EQ.0) THEN
2161 q2s=q1(1+(i-1)*ixq1+(n-1)*nxq1)
2162! J2S=0
2163 ELSEIF(k1.EQ.km1) THEN
2164 q2s=q1(1+(i-1)*ixq1+(km1-1)*kxq1+(n-1)*nxq1)
2165! J2S=0
2166 ELSEIF(k1.EQ.1.OR.k1.EQ.km1-1) THEN
2167 q1a=q1(1+(i-1)*ixq1+(k1-1)*kxq1+(n-1)*nxq1)
2168 q1b=q1(1+(i-1)*ixq1+(k1+0)*kxq1+(n-1)*nxq1)
2169 q2s=ffa(i)*q1a+ffb(i)*q1b
2170! J2S=GGA(I)*Q1A+GGB(I)*Q1B
2171 ELSE
2172 q1a=q1(1+(i-1)*ixq1+(k1-2)*kxq1+(n-1)*nxq1)
2173 q1b=q1(1+(i-1)*ixq1+(k1-1)*kxq1+(n-1)*nxq1)
2174 q1c=q1(1+(i-1)*ixq1+(k1+0)*kxq1+(n-1)*nxq1)
2175 q1d=q1(1+(i-1)*ixq1+(k1+1)*kxq1+(n-1)*nxq1)
2176 q2s=ffa(i)*q1a+ffb(i)*q1b+ffc(i)*q1c+ffd(i)*q1d
2177! J2S=GGA(I)*Q1A+GGB(I)*Q1B+GGC(I)*Q1C+GGD(I)*Q1D
2178 IF(q2s.LT.min(q1b,q1c)) THEN
2179 q2s=min(q1b,q1c)
2180! J2S=0
2181 ELSEIF(q2s.GT.max(q1b,q1c)) THEN
2182 q2s=max(q1b,q1c)
2183! J2S=0
2184 ENDIF
2185 ENDIF
2186 q2(1+(i-1)*ixq2+(k2-1)*kxq2+(n-1)*nxq2)=q2s
2187! J2(1+(I-1)*IXQ2+(K2-1)*KXQ2+(N-1)*NXQ2)=J2S
2188 ENDDO
2189 ENDDO
2190 ENDDO
2191!$OMP END PARALLEL DO
2192
2193 END SUBROUTINE terp3
2194
2251 SUBROUTINE rsearch(IM,KM1,IXZ1,KXZ1,Z1,KM2,IXZ2,KXZ2,Z2,IXL2,KXL2,L2)
2252 IMPLICIT NONE
2253
2254 INTEGER,INTENT(IN) :: IM,KM1,IXZ1,KXZ1,KM2,IXZ2,KXZ2,IXL2,KXL2
2255 INTEGER,INTENT(OUT) :: L2(1+(IM-1)*IXL2+(KM2-1)*KXL2)
2256
2257 REAL(ESMF_KIND_R8),INTENT(IN) :: Z1(1+(IM-1)*IXZ1+(KM1-1)*KXZ1)
2258 REAL(ESMF_KIND_R8),INTENT(IN) :: Z2(1+(IM-1)*IXZ2+(KM2-1)*KXZ2)
2259
2260 INTEGER :: I,K2,L
2261
2262 REAL(ESMF_KIND_R8) :: Z
2263
2264
2265! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2266! FIND THE SURROUNDING INPUT INTERVAL FOR EACH OUTPUT POINT.
2267 DO i=1,im
2268 IF (z1(1+(i-1)*ixz1).LE.z1(1+(i-1)*ixz1+(km1-1)*kxz1)) THEN
2269! INPUT COORDINATE IS MONOTONICALLY ASCENDING.
2270 DO k2=1,km2
2271 z=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
2272 l=0
2273 DO
2274 IF(z.LT.z1(1+(i-1)*ixz1+l*kxz1)) EXIT
2275 l=l+1
2276 IF(l.EQ.km1) EXIT
2277 ENDDO
2278 l2(1+(i-1)*ixl2+(k2-1)*kxl2)=l
2279 ENDDO
2280 ELSE
2281! INPUT COORDINATE IS MONOTONICALLY DESCENDING.
2282 DO k2=1,km2
2283 z=z2(1+(i-1)*ixz2+(k2-1)*kxz2)
2284 l=0
2285 DO
2286 IF(z.GT.z1(1+(i-1)*ixz1+l*kxz1)) EXIT
2287 l=l+1
2288 IF(l.EQ.km1) EXIT
2289 ENDDO
2290 l2(1+(i-1)*ixl2+(k2-1)*kxl2)=l
2291 ENDDO
2292 ENDIF
2293 ENDDO
2294
2295 END SUBROUTINE rsearch
2296
2299 subroutine compute_zh
2300
2301 implicit none
2302
2303 integer :: i,ii, j,k, rc, clb(2), cub(2)
2304
2305 real(esmf_kind_r8), allocatable :: pe0(:), pn0(:)
2306 real(esmf_kind_r8), pointer :: psptr(:,:)
2307 real(esmf_kind_r8), pointer :: zhsfcptr(:,:)
2308 real(esmf_kind_r8), pointer :: zhptr(:,:,:)
2309 real(esmf_kind_r8), pointer :: tptr(:,:,:)
2310 real(esmf_kind_r8), pointer :: qptr(:,:,:)
2311 real(esmf_kind_r8) :: ak, bk, zvir, grd
2312 real(esmf_kind_r8), parameter :: grav = 9.80665
2313 real(esmf_kind_r8), parameter :: rdgas = 287.05
2314 real(esmf_kind_r8), parameter :: rvgas = 461.50
2315
2316 print*,"- COMPUTE HEIGHT"
2317
2318 print*,"- CALL FieldGet FOR SURFACE PRESSURE"
2319 call esmf_fieldget(ps_target_grid, &
2320 computationallbound=clb, &
2321 computationalubound=cub, &
2322 farrayptr=psptr, rc=rc)
2323 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2324 call error_handler("IN FieldGet", rc)
2325
2326 print*,"- CALL FieldGet FOR TERRAIN HEIGHT"
2327 call esmf_fieldget(terrain_target_grid, &
2328 farrayptr=zhsfcptr, rc=rc)
2329 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2330 call error_handler("IN FieldGet", rc)
2331
2332 print*,"- CALL FieldGet FOR HEIGHT"
2333 call esmf_fieldget(zh_target_grid, &
2334 farrayptr=zhptr, rc=rc)
2335 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2336 call error_handler("IN FieldGet", rc)
2337
2338 print*,"- CALL FieldGet FOR TEMPERATURE"
2339 call esmf_fieldget(temp_target_grid, &
2340 farrayptr=tptr, rc=rc)
2341 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2342 call error_handler("IN FieldGet", rc)
2343
2344 do ii = 1, num_tracers
2345 if (trim(tracers(ii)) == "sphum") exit
2346 enddo
2347
2348 print*,"- CALL FieldGet FOR SPECIFIC HUMIDITY"
2349 call esmf_fieldget(tracers_target_grid(ii), &
2350 farrayptr=qptr, rc=rc)
2351 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
2352 call error_handler("IN FieldGet", rc)
2353
2354 grd = grav/rdgas
2355 zvir = rvgas/rdgas - 1.0_esmf_kind_r8
2356
2357 allocate(pe0(levp1_target))
2358 allocate(pn0(levp1_target))
2359
2360 do j = clb(2), cub(2)
2361 do i = clb(1), cub(1)
2362
2363 do k = 1, levp1_target
2364 ak = vcoord_target(k,1)
2365 ak = max(ak, 1.e-9)
2366 bk = vcoord_target(k,2)
2367
2368 pe0(k) = ak + bk*psptr(i,j)
2369 pn0(k) = log(pe0(k))
2370 enddo
2371
2372 zhptr(i,j,1) = zhsfcptr(i,j)
2373
2374 do k = 2, levp1_target
2375 zhptr(i,j,k) = zhptr(i,j,k-1)+tptr(i,j,k-1)*(1.+zvir*qptr(i,j,k-1))* &
2376 (pn0(k-1)-pn0(k))/grd
2377 enddo
2378
2379 enddo
2380 enddo
2381
2382 deallocate(pe0, pn0)
2383
2384 end subroutine compute_zh
2385
2390
2391 implicit none
2392
2393 integer :: i, rc
2394
2395 print*,"- DESTROY TARGET GRID ATMOSPHERIC BEFORE ADJUSTMENT FIELDS."
2396
2397 call esmf_fielddestroy(xwind_b4adj_target_grid, rc=rc)
2398 call esmf_fielddestroy(ywind_b4adj_target_grid, rc=rc)
2399 call esmf_fielddestroy(zwind_b4adj_target_grid, rc=rc)
2400 call esmf_fielddestroy(dzdt_b4adj_target_grid, rc=rc)
2401 call esmf_fielddestroy(ps_b4adj_target_grid, rc=rc)
2402 call esmf_fielddestroy(pres_b4adj_target_grid, rc=rc)
2403 call esmf_fielddestroy(temp_b4adj_target_grid, rc=rc)
2404 call esmf_fielddestroy(terrain_interp_to_target_grid, rc=rc)
2405
2406 do i = 1, num_tracers_input
2407 call esmf_fielddestroy(tracers_b4adj_target_grid(i), rc=rc)
2408 enddo
2409
2410 deallocate(tracers_b4adj_target_grid)
2411
2412 end subroutine cleanup_target_atm_b4adj_data
2413
2417
2419
2420 implicit none
2421
2422 integer :: rc
2423
2424 print*,"- DESTROY LOCAL TARGET GRID ATMOSPHERIC FIELDS."
2425
2426 call esmf_fielddestroy(pres_target_grid, rc=rc)
2427 call esmf_fielddestroy(xwind_target_grid, rc=rc)
2428 call esmf_fielddestroy(ywind_target_grid, rc=rc)
2429 call esmf_fielddestroy(zwind_target_grid, rc=rc)
2430 call esmf_fielddestroy(xwind_s_target_grid, rc=rc)
2431 call esmf_fielddestroy(ywind_s_target_grid, rc=rc)
2432 call esmf_fielddestroy(zwind_s_target_grid, rc=rc)
2433 call esmf_fielddestroy(xwind_w_target_grid, rc=rc)
2434 call esmf_fielddestroy(ywind_w_target_grid, rc=rc)
2435 call esmf_fielddestroy(zwind_w_target_grid, rc=rc)
2436
2438
2439 end subroutine cleanup_all_target_atm_data
2440
2441 end module atmosphere
Read atmospheric data on the input grid.
type(esmf_field), public temp_input_grid
temperature
type(esmf_field), public zwind_input_grid
z-component wind
type(esmf_field), public terrain_input_grid
terrain height
type(esmf_field), public xwind_input_grid
x-component wind
type(esmf_field), public dzdt_input_grid
vert velocity
type(esmf_field), public pres_input_grid
3-d pressure
type(esmf_field), dimension(:), allocatable, public tracers_input_grid
tracers
type(esmf_field), public ps_input_grid
surface pressure
subroutine, public cleanup_input_atm_data
Free up memory associated with atm data.
integer, public lev_input
number of atmospheric layers
type(esmf_field), public ywind_input_grid
y-component wind
integer, public levp1_input
number of atmos layer interfaces
subroutine, public read_input_atm_data(localpet)
Read input grid atmospheric data driver.
Module to hold variables and ESMF fields associated with the target grid atmospheric data.
subroutine, public cleanup_atmosphere_target_data
Free up memory for fields and variables in this module.
type(esmf_field), public u_s_target_grid
U-wind, 'south' edge of grid cell.
type(esmf_field), public ps_target_grid
Surface pressure.
type(esmf_field), public v_s_target_grid
V-wind, 'south' edge of grid cell.
type(esmf_field), public dzdt_target_grid
Vertical velocity.
integer, public levp1_target
Number of vertical levels plus 1.
type(esmf_field), public zh_target_grid
3-d height.
integer, public nvcoord_target
Number of vertical coordinate variables.
type(esmf_field), public u_w_target_grid
U-wind, 'west' edge of grid cell.
type(esmf_field), public temp_target_grid
Temperautre.
type(esmf_field), dimension(:), allocatable, public tracers_target_grid
Tracers.
integer, public lev_target
Number of vertical levels.
type(esmf_field), public qnwfa_climo_target_grid
Number concentration of water friendly aerosols.
type(esmf_field), public qnifa_climo_target_grid
Number concentration of ice friendly aerosols.
type(esmf_field), public delp_target_grid
Pressure thickness.
type(esmf_field), public v_w_target_grid
V-wind, 'west' edge of grid cell.
real(esmf_kind_r8), dimension(:,:), allocatable, public vcoord_target
Vertical coordinate.
Process atmospheric fields.
type(esmf_field) zwind_target_grid
z-component wind, grid box center
subroutine create_atm_esmf_fields
Create target grid field objects.
type(esmf_field) ywind_s_target_grid
y-component wind, 'south' edge
type(esmf_field) zwind_s_target_grid
z-component wind, 'south' edge
subroutine vintg_wam(year, month, day, hour, pf)
Vertically extend model top into thermosphere for whole atmosphere model.
type(esmf_field) ywind_w_target_grid
y-component wind, 'west' edge
type(esmf_field) ywind_b4adj_target_grid
y-component wind, before vert adj
type(esmf_field) xwind_b4adj_target_grid
x-component wind, before vert adj
subroutine horiz_interp_thomp_mp_climo
Horizontally interpolate thompson microphysics data to the target model grid.
type(esmf_field) pres_target_grid
3-d pressure
type(esmf_field) pres_b4adj_target_grid
3-d pres before terrain adj
subroutine newpr1(localpet)
Computes 3-D pressure given an adjusted surface pressure.
subroutine vintg_thomp_mp_climo
Vertically interpolate atmospheric fields to target FV3 grid.
type(esmf_field) temp_b4adj_target_grid
temp before vert adj
subroutine, public read_vcoord_info
Reads model vertical coordinate definition file (as specified by namelist variable vcoord_file_target...
subroutine newps(localpet)
Computes adjusted surface pressure given a new terrain height.
subroutine cleanup_target_atm_b4adj_data
Cleanup atmospheric field (before adjustment) objects.
type(esmf_field) terrain_interp_to_target_grid
Input grid terrain interpolated to target grid.
subroutine rsearch(im, km1, ixz1, kxz1, z1, km2, ixz2, kxz2, z2, ixl2, kxl2, l2)
Search for a surrounding real interval.
type(esmf_field) xwind_s_target_grid
x-component wind, 'south' edge
type(esmf_field) ywind_target_grid
y-component wind, grid box center
subroutine create_atm_b4adj_esmf_fields
Create target grid field objects to hold data before vertical interpolation.
type(esmf_field) ps_b4adj_target_grid
sfc pres before terrain adj
subroutine compute_zh
Compute vertical level height.
subroutine, public atmosphere_driver(localpet)
Driver routine to process for atmospheric fields.
subroutine terp3(im, ixz1, ixq1, ixz2, ixq2, nm, nxq1, nxq2, km1, kxz1, kxq1, z1, q1, km2, kxz2, kxq2, z2, q2)
Cubically interpolate in one dimension.
type(esmf_field) xwind_w_target_grid
x-component wind, 'west' edge
type(esmf_field) dzdt_b4adj_target_grid
vertical vel before vert adj
type(esmf_field) xwind_target_grid
x-component wind, grid box center
type(esmf_field) thomp_pres_climo_b4adj_target_grid
pressure of each level on target grid
type(esmf_field) zwind_b4adj_target_grid
z-component wind, before vert adj
subroutine vintg
Vertically interpolate upper-air fields.
type(esmf_field) zwind_w_target_grid
z-component wind, 'west' edge
type(esmf_field), dimension(:), allocatable tracers_b4adj_target_grid
tracers before vert adj
subroutine convert_winds_to_uv
Convert 3-d component winds to u and v.
type(esmf_field) qnifa_climo_b4adj_target_grid
number concentration of ice friendly aerosols before vert adj
type(esmf_field) qnwfa_climo_b4adj_target_grid
number concentration of water friendly aerosols before vert adj
subroutine cleanup_all_target_atm_data
Cleanup target grid atmospheric field objects.
Sets up the ESMF grid objects for the input data grid and target FV3 grid.
Definition model_grid.F90:9
type(esmf_field), public latitude_w_target_grid
latitude of 'west' edge of grid box, target grid
type(esmf_field), public longitude_w_target_grid
longitude of 'west' edge of grid box, target grid
type(esmf_field), public terrain_target_grid
terrain height target grid
type(esmf_grid), public target_grid
target grid esmf grid object.
type(esmf_field), public longitude_s_target_grid
longitude of 'south' edge of grid box, target grid
type(esmf_field), public latitude_s_target_grid
latitude of 'south' edge of grid box, target grid
This module contains code to read the setup namelist file, handle the varmap file for GRIB2 data,...
integer, public num_tracers_input
Number of atmospheric tracers in input file.
integer, public regional
For regional target grids.
character(len=500), public wam_parm_file
Full path to msis21.parm for WAM initialization.
character(len=20), dimension(max_tracers), public tracers
Name of each atmos tracer to be processed.
integer, public num_tracers
Number of atmospheric tracers to be processed.
logical, public use_thomp_mp_climo
When true, read and process Thompson MP climatological tracers.
logical, public wam_cold_start
When true, cold start for whole atmosphere model.
integer, public cycle_mon
Cycle month.
integer, public cycle_day
Cycle day.
integer, public cycle_hour
Cycle hour.
character(len=500), public atm_weight_file
File containing pre-computed weights to horizontally interpolate atmospheric fields.
character(len=500), public vcoord_file_target_grid
Vertical coordinate definition file.
integer, public cycle_year
Cycle year.
Module to read the Thompson climatological MP data file and set up the associated esmf field and grid...
type(esmf_field), public qnifa_climo_input_grid
number concentration of ice friendly nuclei.
subroutine, public cleanup_thomp_mp_climo_input_data
Free up memory associated with this module.
type(esmf_field), public qnwfa_climo_input_grid
number concentration of water friendly nuclei.
type(esmf_field), public thomp_pres_climo_input_grid
3-d pressure of the Thompson climo data points
subroutine, public read_thomp_mp_climo_data
Read Thompson climatological MP data file and time interpolate data to current cycle time.
integer, public lev_thomp_mp_climo
number of vert lvls of Thompson climo data
subroutine gettemp(iday, xlat, pr, np, pf, temp, n_o, n_o2, n_n2)
Routine that computes temperature and neutral density values utilizing MSIS 2.1.