forked from wrf-model/WRF
-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathmodule_initialize_real.F
9187 lines (7767 loc) · 411 KB
/
module_initialize_real.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!REAL:MODEL_LAYER:INITIALIZATION
#ifndef VERT_UNIT
! This MODULE holds the routines which are used to perform various initializations
! for the individual domains, specifically for the Eulerian, mass-based coordinate.
!-----------------------------------------------------------------------
MODULE module_initialize_real
USE module_bc
USE module_configure
USE module_domain
USE module_io_domain
USE module_model_constants
USE module_state_description
USE module_timing
USE module_soil_pre
USE module_date_time
USE module_llxy
USE module_polarfft
#ifdef DM_PARALLEL
USE module_dm
USE module_comm_dm, ONLY : &
HALO_EM_INIT_1_sub &
,HALO_EM_INIT_2_sub &
,HALO_EM_INIT_3_sub &
,HALO_EM_INIT_4_sub &
,HALO_EM_INIT_5_sub &
,HALO_EM_INIT_6_sub &
,HALO_EM_VINTERP_UV_1_sub
#endif
REAL , SAVE :: p_top_save
INTEGER :: internal_time_loop
CONTAINS
!-------------------------------------------------------------------
SUBROUTINE init_domain ( grid )
IMPLICIT NONE
! Input space and data. No gridded meteorological data has been stored, though.
! TYPE (domain), POINTER :: grid
TYPE (domain) :: grid
! Local data.
INTEGER :: idum1, idum2
CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
CALL init_domain_rk( grid &
!
#include "actual_new_args.inc"
!
)
END SUBROUTINE init_domain
!-------------------------------------------------------------------
SUBROUTINE init_domain_rk ( grid &
!
#include "dummy_new_args.inc"
!
)
USE module_optional_input
USE module_radiation_driver, ONLY: cal_cldfra3
USE module_dm, ONLY : wrf_dm_max_real
use module_madwrf, only : Init_madwrf_clouds, Init_madwrf_tracers
IMPLICIT NONE
! Input space and data. No gridded meteorological data has been stored, though.
! TYPE (domain), POINTER :: grid
TYPE (domain) :: grid
#include "dummy_new_decl.inc"
TYPE (grid_config_rec_type) :: config_flags
! Local domain indices and counters.
INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
INTEGER :: loop , num_seaice_changes
INTEGER :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte, &
ips, ipe, jps, jpe, kps, kpe, &
i, j, k, kk
INTEGER :: imsx, imex, jmsx, jmex, kmsx, kmex, &
ipsx, ipex, jpsx, jpex, kpsx, kpex, &
imsy, imey, jmsy, jmey, kmsy, kmey, &
ipsy, ipey, jpsy, jpey, kpsy, kpey
INTEGER :: ns
! Local data
INTEGER :: error
INTEGER :: im, num_3d_m, num_3d_s
REAL :: B1, B2, B3, B4, B5
REAL :: p_surf, p_level
REAL :: cof1, cof2
REAL :: qvf , qvf1 , qvf2 , qtot, pd_surf
REAL :: p00 , t00 , a , tiso, p_strat, a_strat
REAL :: hold_znw , ptemp
REAL :: vap_pres_mb , sat_vap_pres_mb
LOGICAL :: were_bad
LOGICAL :: stretch_grid, dry_sounding, debug
INTEGER :: IICOUNT, icount
REAL :: p_top_requested , temp
INTEGER :: num_metgrid_levels
INTEGER, PARAMETER :: num_wif_levels_default = 30
INTEGER :: aer_init_opt
REAL , DIMENSION(max_eta) :: eta_levels
INTEGER :: auto_levels_opt
REAL :: max_dz, dzbot, dzstretch_s, dzstretch_u, z1
INTEGER:: i_end, j_end
REAL, ALLOCATABLE, DIMENSION(:):: temp_P, temp_T, temp_R, temp_Dz
REAL, ALLOCATABLE, DIMENSION(:):: temp_Qv, temp_Qc, temp_Qi, temp_Qs
REAL, ALLOCATABLE, DIMENSION(:):: temp_CF, temp_Nc, temp_Ni
REAL:: max_relh, temp_xland, gridkm
LOGICAL :: debug_flag = .FALSE.
REAL :: dclat
! INTEGER , PARAMETER :: nl_max = 1000
! REAL , DIMENSION(nl_max) :: grid%dn
integer::oops1,oops2
REAL :: zap_close_levels
INTEGER :: force_sfc_in_vinterp
INTEGER :: interp_type , lagrange_order , extrap_type , t_extrap_type
INTEGER :: linear_interp
LOGICAL :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
LOGICAL :: we_have_tavgsfc , we_have_tsk
INTEGER :: lev500 , loop_count
REAL :: zl , zu , pl , pu , z500 , dz500 , tvsfc , dpmu
REAL :: pfu, pfd, phm
LOGICAL , PARAMETER :: want_full_levels = .TRUE.
LOGICAL , PARAMETER :: want_half_levels = .FALSE.
CHARACTER (LEN=256) :: a_message, mminlu
REAL :: max_mf
! Excluded middle.
LOGICAL :: any_valid_points
INTEGER :: i_valid , j_valid
! Vert interpolation in WRF
INTEGER :: k_max_p , k_min_p
!-- Carsel and Parrish [1988]
REAL , DIMENSION(100) :: lqmi
REAL , DIMENSION(100) :: thickness , levels
REAL :: t_start , t_end
REAL , ALLOCATABLE , DIMENSION(:,:) :: clat_glob
! added for multiple specified sets of eta_levels with vertical grid nesting
INTEGER :: ks, ke, id
LOGICAL :: vnest !T if using vertical nesting with vet_refine_method=2, otherwise F
INTEGER :: j_save
INTEGER :: change_soil, change_soilw, iforce
REAL:: temp_rho
LOGICAL :: wif_upside_down = .FALSE.
! Test on consistency between namelist settings and the available data from geogrid.
INTEGER :: geogrid_flag_error
! Vertical pressure checks
REAL :: press_above, press_below
! Dimension information stored in grid data structure.
CALL cpu_time(t_start)
CALL get_ijk_from_grid ( grid , &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
imsx, imex, jmsx, jmex, kmsx, kmex, &
ipsx, ipex, jpsx, jpex, kpsx, kpex, &
imsy, imey, jmsy, jmey, kmsy, kmey, &
ipsy, ipey, jpsy, jpey, kpsy, kpey )
its = ips ; ite = ipe ; jts = jps ; jte = jpe ; kts = kps ; kte = kpe
CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
! Check flags from geogrid and the various model settings for any inconsistency.
! The optional data from geogrid mostly started with the v4.0 release. Older
! WPS data could have the required fields. Assume users know what they are
! doing when they bring old data into the newer real program ...
IF ( grid%v4_metgrid ) THEN
geogrid_flag_error = 0
IF ( ( config_flags%topo_wind .EQ. 1 ) .AND. ( flag_var_sso .EQ. 0 ) ) THEN
CALL wrf_message ( '----- ERROR: topo_wind = 1 AND flag_var_sso = 0 ' )
geogrid_flag_error = geogrid_flag_error + 1
END IF
IF ( ( config_flags%sf_lake_physics .EQ. 1 ) .AND. ( flag_lake_depth .EQ. 0 ) ) THEN
CALL wrf_message ( '----- ERROR: sf_lake_physics = 1 AND flag_lake_depth = 0 ' )
geogrid_flag_error = geogrid_flag_error + 1
END IF
IF ( ( config_flags%sf_surface_physics .EQ. pxlsmscheme ) .AND. ( flag_imperv .EQ. 0 ) ) THEN
CALL wrf_message ( '----- ERROR: sf_surface_physics = 7 AND flag_imperv = 0 ' )
geogrid_flag_error = geogrid_flag_error + 1
END IF
IF ( ( config_flags%sf_surface_physics .EQ. pxlsmscheme ) .AND. ( flag_canfra .EQ. 0 ) ) THEN
CALL wrf_message ( '----- ERROR: sf_surface_physics = 7 AND flag_canfra = 0 ' )
geogrid_flag_error = geogrid_flag_error + 1
END IF
IF ( ( config_flags%mp_physics .EQ. thompsonaero ) .AND. &
( config_flags%dust_emis .EQ. 1 ) .AND. ( flag_erod .EQ. 0 ) ) THEN
CALL wrf_message ( '----- ERROR: mp=28 AND dust_emis= 1 AND flag_erod = 0 ' )
geogrid_flag_error = geogrid_flag_error + 1
END IF
IF ( ( config_flags%mp_physics .EQ. thompsonaero ) .AND. &
( config_flags%dust_emis .EQ. 1 ) .AND. ( flag_clayfrac .EQ. 0 ) ) THEN
CALL wrf_message ( '----- ERROR: mp=28 AND dust_emis= 1 AND flag_clayfrac = 0 ' )
geogrid_flag_error = geogrid_flag_error + 1
END IF
IF ( ( config_flags%mp_physics .EQ. thompsonaero ) .AND. &
( config_flags%dust_emis .EQ. 1 ) .AND. ( flag_sandfrac .EQ. 0 ) ) THEN
CALL wrf_message ( '----- ERROR: mp=28 AND dust_emis= 1 AND flag_sandfrac = 0 ' )
geogrid_flag_error = geogrid_flag_error + 1
END IF
IF ( geogrid_flag_error .GT. 0 ) THEN
CALL wrf_error_fatal ('Either modify the namelist settings, or rebuild the geogrid/metgrid data' )
END IF
! Geogrid flags that are not yet used: FLAG_FRC_URB2D FLAG_LAI12M FLAG_URB_PARAM
END IF
! Would the user prefer to forego the use of the level of max winds, or the
! tropopause level data? This is an option the user may select. While the
! additional data is able to provide good information (such as a better
! resolution of the jet, a better kink for the tropopause), there are
! horizontal gradients that are introduced. Near a boundary, these gradients
! would be permanent (due to their inclusion in the LBC file). To turn "off"
! the use of the max wind/trop data, set the flags for those levels to zero.
flag_pmaxw = 0
flag_pmaxwnn = 0
flag_ptrop = 0
flag_ptropnn = 0
IF ( ( config_flags%use_maxw_level .EQ. 0 ) .AND. &
( ( flag_tmaxw .EQ. 1 ) .OR. ( flag_umaxw .EQ. 1 ) .OR. ( flag_vmaxw .EQ. 1 ) .OR. ( flag_hgtmaxw .EQ. 1 ) ) ) THEN
flag_tmaxw = 0
flag_umaxw = 0
flag_vmaxw = 0
flag_hgtmaxw = 0
CALL wrf_debug ( 0 , 'Turning off use of MAX WIND level data in vertical interpolation' )
END IF
IF ( ( config_flags%use_trop_level .EQ. 0 ) .AND. &
( ( flag_ttrop .EQ. 1 ) .OR. ( flag_utrop .EQ. 1 ) .OR. ( flag_vtrop .EQ. 1 ) .OR. ( flag_hgttrop .EQ. 1 ) ) ) THEN
flag_ttrop = 0
flag_utrop = 0
flag_vtrop = 0
flag_hgttrop = 0
CALL wrf_debug ( 0 , 'Turning off use of TROPOPAUSE level data in vertical interpolation' )
END IF
! Lake Mask and depth assignment
CALL nl_get_iswater ( grid%id , grid%iswater )
CALL nl_get_islake ( grid%id , grid%islake )
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
IF ( grid%lu_index(i,j) .NE. grid%islake ) THEN
grid%lakemask(i,j) = 0
ELSE
grid%lakemask(i,j) = 1
END IF
END DO
END DO
IF ( grid%sf_lake_physics .EQ. 1 ) THEN
grid%lake_depth_flag = flag_lake_depth
IF ( flag_lake_depth .EQ. 0 ) THEN
CALL wrf_message ( " Warning: Please rerun WPS to get lake_depth information for lake model" )
! Set lake depth over the ocean to be -3 m, and set the lake depth over land to be -2 m.
ELSE
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
IF ( ( grid%lu_index(i,j) .NE. grid%islake ) .AND. ( grid%lu_index(i,j) .NE. grid%iswater ) ) THEN
grid%lake_depth(i,j) = -2
ELSE IF ( grid%lu_index(i,j) .NE. grid%islake ) THEN
grid%lake_depth(i,j) = -3
END IF
END DO
END DO
END IF
END IF
grid%bathymetry_flag = flag_bathymetry
IF ( flag_bathymetry .EQ. 0 ) THEN
IF ( grid%shalwater_z0 .EQ. 1 ) THEN
CALL wrf_message ( " Warning: No bathymetry data found for shallow water roughness model." )
IF ( grid%shalwater_depth .LE. 0.0 ) THEN
CALL wrf_message ( " Warning: shalwater_depth must be greater than 0.0 for WRF to run." )
END IF
END IF
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
grid%water_depth(i,j) = -4.0
END DO
END DO
ELSE
CALL wrf_message ( " Bathymetry dataset from GEBCO Compilation Group. Please acknowledge the following in presentations and publications: GEBCO Compilation Group (2021) GEBCO 2021 Grid (doi:10.5285/c6612cbe-50b3-0cff-e053-6c86abc09f8f)." )
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
grid%water_depth(i,j) = grid%bathymetry(i,j)
! Get depth of lake based on height of water surface:
IF ( grid%lu_index(i,j) .EQ. grid%islake ) THEN
grid%water_depth(i,j) = grid%bathymetry(i,j) - grid%ht_gc(i,j)
END IF
! Depth is positive:
grid%water_depth(i,j) = -grid%water_depth(i,j)
! Set land cells to -10
IF ( ( grid%lu_index(i,j) .NE. grid%islake ) .AND. ( grid%lu_index(i,j) .NE. grid%iswater ) ) THEN
grid%water_depth(i,j) = -2.0
ELSE ! Water cells:
! Find any water cells with negative (originally positive) values...
! ... indicative of mis-match of bathymetry and land mask.
IF (grid%water_depth(i,j) .LT. 0.1) THEN
grid%water_depth(i,j) = 0.1
END IF
END IF
END DO
END DO
END IF
! Send out a quick message about the time steps based on the map scale factors.
IF ( ( internal_time_loop .EQ. 1 ) .AND. ( grid%id .EQ. 1 ) .AND. &
( .NOT. config_flags%polar ) ) THEN
max_mf = grid%msft(its,jts)
DO j=jts,MIN(jde-1,jte)
DO i=its,MIN(ide-1,ite)
max_mf = MAX ( max_mf , grid%msft(i,j) )
END DO
END DO
#if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
max_mf = wrf_dm_max_real ( max_mf )
#endif
WRITE ( a_message , FMT='(A,F5.2,A)' ) 'Max map factor in domain 1 = ',max_mf, &
'. Scale the dt in the model accordingly.'
CALL wrf_message ( a_message )
END IF
! Check to see if the boundary conditions are set properly in the namelist file.
! This checks for sufficiency and redundancy.
CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
! Some sort of "this is the first time" initialization. Who knows.
grid%step_number = 0
grid%itimestep=0
! Pull in the info in the namelist to compare it to the input data.
grid%real_data_init_type = model_config_rec%real_data_init_type
! To define the base state, we call a USER MODIFIED routine to set the three
! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K),
! and A (temperature difference, from 1000 mb to 300 mb, K).
CALL const_module_initialize ( p00 , t00 , a , tiso , p_strat , a_strat )
! Save these constants to write out in model output file
grid%t00 = t00
grid%p00 = p00
grid%tlp = a
grid%tiso = tiso
grid%p_strat = p_strat
grid%tlp_strat = a_strat
! Are there any hold-ups to us bypassing the middle of the domain? These
! holdups would be situations where we need data in the middle of the domain.
! FOr example, if this si the first time period, we need the full domain
! processed for ICs. Also, if there is some sort of gridded FDDA turned on, or
! if the SST update is activated, then we can't just blow off the middle of the
! domain all willy-nilly. Other cases of these hold-ups? Sure - what if the
! user wants to smooth the CG topo, we need several rows and columns available.
! What if the lat/lon proj is used, then we need to run a spectral filter on
! the topo. Both are killers when trying to ignore data in the middle of the
! domain.
! If hold_ups = .F., then there are no hold-ups to excluding the middle
! domain processing. If hold_ups = .T., then there are hold-ups, and we
! must process the middle of the domain.
hold_ups = ( internal_time_loop .EQ. 1 ) .OR. &
( config_flags%grid_fdda .NE. 0 ) .OR. &
( config_flags%sst_update .EQ. 1 ) .OR. &
( config_flags%qna_update .EQ. 1 ) .OR. &
( config_flags%all_ic_times ) .OR. &
( config_flags%polar )
! There are a few checks that we need to do when the input data comes in with the middle
! excluded by WPS.
IF ( flag_excluded_middle .NE. 0 ) THEN
! If this time period of data from WPS has the middle excluded, it had better be OK for
! us to have a hole.
IF ( hold_ups ) THEN
WRITE ( a_message,* ) 'None of the following are allowed to be TRUE : '
CALL wrf_message ( a_message )
WRITE ( a_message,* ) ' ( internal_time_loop .EQ. 1 ) ', ( internal_time_loop .EQ. 1 )
CALL wrf_message ( a_message )
WRITE ( a_message,* ) ' ( config_flags%grid_fdda .NE. 0 ) ', ( config_flags%grid_fdda .NE. 0 )
CALL wrf_message ( a_message )
WRITE ( a_message,* ) ' ( config_flags%sst_update .EQ. 1 ) ', ( config_flags%sst_update .EQ. 1 )
CALL wrf_message ( a_message )
WRITE ( a_message,* ) ' ( config_flags%qna_update .EQ. 1 ) ', ( config_flags%qna_update .EQ. 1 )
CALL wrf_message ( a_message )
WRITE ( a_message,* ) ' ( config_flags%all_ic_times ) ', ( config_flags%all_ic_times )
CALL wrf_message ( a_message )
WRITE ( a_message,* ) ' ( config_flags%smooth_cg_topo ) ', ( config_flags%smooth_cg_topo )
CALL wrf_message ( a_message )
WRITE ( a_message,* ) ' ( config_flags%polar ) ', ( config_flags%polar )
CALL wrf_message ( a_message )
WRITE ( a_message,* ) 'Problems, we cannot have excluded middle data from WPS'
CALL wrf_error_fatal ( a_message )
END IF
! Make sure that the excluded middle data from metgrid is "wide enough". We only have to check
! when the excluded middle was actually used in WPS.
IF ( config_flags%spec_bdy_width .GT. flag_excluded_middle ) THEN
WRITE ( a_message,* ) 'The WRF &bdy_control namelist.input spec_bdy_width = ', config_flags%spec_bdy_width
CALL wrf_message ( a_message )
WRITE ( a_message,* ) 'The WPS &metgrid namelist.wps process_only_bdy width = ',flag_excluded_middle
CALL wrf_message ( a_message )
WRITE ( a_message,* ) 'WPS process_only_bdy must be >= WRF spec_bdy_width'
CALL wrf_error_fatal ( a_message )
END IF
END IF
em_width = config_flags%spec_bdy_width
! We need to find if there are any valid non-excluded-middle points in this
! tile. If so, then we need to hang on to a valid i,j location.
any_valid_points = .false.
find_valid : DO j = jts,jte
DO i = its,ite
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
any_valid_points = .true.
i_valid = i
j_valid = j
EXIT find_valid
END DO
END DO find_valid
! Replace traditional seaice field with optional seaice (AFWA source)
IF ( flag_icefrac .EQ. 1 ) THEN
DO j=jts,MIN(jde-1,jte)
DO i=its,MIN(ide-1,ite)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%xice(i,j) = grid%icefrac_gc(i,j)
END DO
END DO
END IF
! Replace traditional seaice field with optional seaice percent (AFWA source)
IF ( flag_icepct .EQ. 1 ) THEN
DO j=jts,MIN(jde-1,jte)
DO i=its,MIN(ide-1,ite)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%xice(i,j) = grid%icepct(i,j)/100.
END DO
END DO
END IF
! Fix the snow (water equivalent depth, kg/m^2) and the snowh (physical snow
! depth, m) fields.
IF ( ( flag_snow .EQ. 0 ) .AND. ( flag_snowh .EQ. 0 ) ) THEN
DO j=jts,MIN(jde-1,jte)
DO i=its,MIN(ide-1,ite)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%snow(i,j) = 0.
grid%snowh(i,j) = 0.
END DO
END DO
ELSE IF ( ( flag_snow .EQ. 0 ) .AND. ( flag_snowh .EQ. 1 ) ) THEN
DO j=jts,MIN(jde-1,jte)
DO i=its,MIN(ide-1,ite)
! ( m -> kg/m^2 ) & ( reduce to liquid, 5:1 ratio )
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%snow(i,j) = grid%snowh(i,j) * 1000. / 5.
END DO
END DO
ELSE IF ( ( flag_snow .EQ. 1 ) .AND. ( flag_snowh .EQ. 0 ) ) THEN
DO j=jts,MIN(jde-1,jte)
DO i=its,MIN(ide-1,ite)
! ( kg/m^2 -> m) & ( liquid to snow depth, 5:1 ratio )
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%snowh(i,j) = grid%snow(i,j) / 1000. * 5.
END DO
END DO
END IF
! For backward compatibility, we might need to assign the map factors from
! what they were, to what they are.
IF ( ( config_flags%polar ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
DO j=max(jds+1,jts),min(jde-1,jte)
DO i=its,min(ide-1,ite)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
END DO
END DO
IF(jts == jds) THEN
DO i=its,ite
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%msfvx(i,jts) = 0.
grid%msfvx_inv(i,jts) = 0.
END DO
END IF
IF(jte == jde) THEN
DO i=its,ite
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%msfvx(i,jte) = 0.
grid%msfvx_inv(i,jte) = 0.
END DO
END IF
ELSE IF ( ( .NOT. config_flags%polar ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
IF ( grid%msfvx(its,jts) .EQ. 0 ) THEN
CALL wrf_error_fatal ( 'Maybe this is a global domain, but the polar flag was not set in the bdy_control namelist.' )
END IF
DO j=jts,min(jde,jte)
DO i=its,min(ide-1,ite)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
END DO
END DO
ELSE IF ( ( config_flags%polar ) .AND. ( flag_mf_xy .NE. 1 ) ) THEN
CALL wrf_error_fatal ( 'Older metgrid data cannot initialize a global domain' )
ENDIF
! Check to see what available surface temperatures we have.
IF ( flag_tavgsfc .EQ. 1 ) THEN
we_have_tavgsfc = .TRUE.
ELSE
we_have_tavgsfc = .FALSE.
END IF
IF ( flag_tsk .EQ. 1 ) THEN
we_have_tsk = .TRUE.
ELSE
we_have_tsk = .FALSE.
END IF
IF ( config_flags%use_tavg_for_tsk ) THEN
IF ( we_have_tsk .OR. we_have_tavgsfc ) THEN
! we are OK
ELSE
CALL wrf_error_fatal ( 'We either need TSK or TAVGSFC, verify these fields are coming from WPS' )
END IF
! Since we require a skin temperature in the model, we can use the average 2-m temperature if provided.
IF ( we_have_tavgsfc ) THEN
DO j=jts,min(jde,jte)
DO i=its,min(ide-1,ite)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%tsk(i,j) = grid%tavgsfc(i,j)
END DO
END DO
END IF
END IF
IF (config_flags%slucm_distributed_drag) THEN
CALL wrf_message('Adding zero-plane displacement height to topography')
DO j = jts, MIN(jde - 1, jte)
DO i = its, MIN(ide - 1, ite)
IF (grid%zd_urb2d(i, j) > 0) grid%ht_gc(i, j) = grid%ht_gc(i, j) + grid%zd_urb2d(i, j)
END DO
END DO
END IF
! Is there any vertical interpolation to do? The "old" data comes in on the correct
! vertical locations already.
IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
num_metgrid_levels = grid%num_metgrid_levels
IF ( config_flags%nest_interp_coord .EQ. 1 ) THEN
! At the location of maximum pressure in the column, get the temperature and height. These
! will be written out and could be used for vertical interpolation - to avoid extrapolation.
! Hey, we can also do minimum values, too.
DO j=jts,jte
DO i=its,ite
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%max_p(i,j) = grid%p_gc(i,1,j)
k_max_p = 1
IF ( grid%p_gc(i,2,j) .GT. grid%max_p(i,j) ) THEN
grid%max_p(i,j) = grid%p_gc(i,2,j)
k_max_p = 2
ELSE IF ( grid%p_gc(i,num_metgrid_levels,j) .GT. grid%max_p(i,j) ) THEN
grid%max_p(i,j) = grid%p_gc(i,num_metgrid_levels,j)
k_max_p = num_metgrid_levels
END IF
grid%t_max_p(i,j) = grid%t_gc(i,k_max_p,j)
grid%ght_max_p(i,j) = grid%ght_gc(i,k_max_p,j)
grid%min_p(i,j) = grid%p_gc(i,num_metgrid_levels,j)
k_min_p = num_metgrid_levels
IF ( grid%p_gc(i,2,j) .LT. grid%min_p(i,j) ) THEN
grid%min_p(i,j) = grid%p_gc(i,2,j)
k_min_p = 2
END IF
grid%t_min_p(i,j) = grid%t_gc(i,k_min_p,j)
grid%ght_min_p(i,j) = grid%ght_gc(i,k_min_p,j)
END DO
END DO
END IF
! If this is data from the PINTERP program, it is emulating METGRID output.
! One of the caveats of this data is the way that the vertical structure is
! handled. We take the k=1 level and toss it (it is disposable), and we
! swap in the surface data. This is done for all of the 3d fields about
! which we show some interest: u, v, t, rh, ght, and p. For u, v, and rh,
! we assume no interesting vertical structure, and just assign the 1000 mb
! data. We directly use the 2-m temp for surface temp. We use the surface
! pressure field and the topography elevation for the lowest level of
! pressure and height, respectively.
IF ( flag_pinterp .EQ. 1 ) THEN
WRITE ( a_message , * ) 'Data from P_INTERP program, filling k=1 level with artificial surface fields.'
CALL wrf_message ( a_message )
DO j=jts,jte
DO i=its,ite
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%u_gc(i,1,j) = grid%u_gc(i,2,j)
grid%v_gc(i,1,j) = grid%v_gc(i,2,j)
grid%rh_gc(i,1,j) = grid%rh_gc(i,2,j)
grid%t_gc(i,1,j) = grid%t2(i,j)
grid%ght_gc(i,1,j) = grid%ht(i,j)
grid%p_gc(i,1,j) = grid%psfc(i,j)
END DO
END DO
flag_psfc = 0
END IF
! Variables that are named differently between SI and WPS.
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%tsk(i,j) = grid%tsk_gc(i,j)
grid%tmn(i,j) = grid%tmn_gc(i,j)
grid%xlat(i,j) = grid%xlat_gc(i,j)
grid%xlong(i,j) = grid%xlong_gc(i,j)
grid%ht(i,j) = grid%ht_gc(i,j)
END DO
END DO
! A user could request that the most coarse grid has the
! topography along the outer boundary smoothed. This smoothing
! is similar to the coarse/nest interface. The outer rows and
! cols come from the existing large scale topo, and then the
! next several rows/cols are a linear ramp of the large scale
! model and the hi-res topo from WPS. We only do this for the
! coarse grid since we are going to make the interface consistent
! in the model betwixt the CG and FG domains.
! An important point is to inform the user if their request cannot
! be satisfied. Do not skip over this quietly.
IF ( ( config_flags%smooth_cg_topo ) .AND. &
( internal_time_loop .EQ. 1 ) .AND. &
( grid%id .EQ. 1 ) .AND. &
( flag_soilhgt .NE. 1) ) THEN
CALL wrf_message (' --- ERROR: NML option smooth_cg_topo=T')
CALL wrf_message (' But found no soil elevation / terrain / topography data in metgrid files')
CALL wrf_message (' The field SOILHGT is required when smoothing the CG topography on d01')
CALL wrf_error_fatal(' If using ERA5 data, possibly need to add more time invariant fields')
END IF
IF ( ( config_flags%smooth_cg_topo ) .AND. &
( internal_time_loop .EQ. 1 ) .AND. &
( grid%id .EQ. 1 ) .AND. &
( flag_soilhgt .EQ. 1) ) THEN
CALL blend_terrain ( grid%toposoil , grid%ht , &
ids , ide , jds , jde , 1 , 1 , &
ims , ime , jms , jme , 1 , 1 , &
ips , ipe , jps , jpe , 1 , 1 )
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
grid%ht_smooth(i,j) = grid%ht(i,j)
END DO
END DO
ELSE IF ( ( config_flags%smooth_cg_topo ) .AND. &
( internal_time_loop .NE. 1 ) .AND. &
( grid%id .EQ. 1 ) .AND. &
( flag_soilhgt .EQ. 1) ) THEN
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
grid%ht(i,j) = grid%ht_smooth(i,j)
END DO
END DO
END IF
! Filter the input topography if this is a global domain.
IF ( ( config_flags%polar ) .AND. ( grid%fft_filter_lat .GT. 90 ) ) THEN
CALL wrf_error_fatal ( 'If the polar boundary condition is used, then fft_filter_lat must be set in namelist.input' )
END IF
IF ( ( config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( config_flags%polar ) ) THEN
#if 1
dclat = 90./REAL(jde-jds) !(0.5 * 180/ny)
DO j = jts, MIN(jte,jde-1)
DO k = kts, kte
DO i = its, MIN(ite,ide-1)
grid%t_2(i,k,j) = 1.
END DO
END DO
DO i = its, MIN(ite,ide-1)
grid%t_2(i,1,j) = grid%ht(i,j)
grid%sr(i,j) = grid%ht(i,j)
END DO
END DO
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! WARNING: this might present scaling issues on very large numbers of processors
ALLOCATE( clat_glob(ids:ide,jds:jde) )
CALL wrf_patch_to_global_real ( grid%clat, clat_glob, grid%domdesc, 'xy', 'xy', &
ids, ide, jds, jde, 1, 1, &
ims, ime, jms, jme, 1, 1, &
its, ite, jts, jte, 1, 1 )
CALL wrf_dm_bcast_real ( clat_glob , (ide-ids+1)*(jde-jds+1) )
grid%clat_xxx(ipsx:ipex,jpsx:jpex) = clat_glob(ipsx:ipex,jpsx:jpex)
find_j_index_of_fft_filter : DO j = jds , jde-1
IF ( ABS(clat_glob(ids,j)) .LE. config_flags%fft_filter_lat ) THEN
j_save = j
EXIT find_j_index_of_fft_filter
END IF
END DO find_j_index_of_fft_filter
CALL wrf_patch_to_global_real ( grid%msft, clat_glob, grid%domdesc, 'xy', 'xy', &
ids, ide, jds, jde, 1, 1, &
ims, ime, jms, jme, 1, 1, &
its, ite, jts, jte, 1, 1 )
CALL wrf_dm_bcast_real ( clat_glob , (ide-ids+1)*(jde-jds+1) )
grid%mf_fft = clat_glob(ids,j_save)
grid%mf_xxx(ipsx:ipex,jpsx:jpex) = clat_glob(ipsx:ipex,jpsx:jpex)
DEALLOCATE( clat_glob )
#else
find_j_index_of_fft_filter : DO j = jds , jde-1
IF ( ABS(grid%clat(ids,j)) .LE. config_flags%fft_filter_lat ) THEN
j_save = j
EXIT find_j_index_of_fft_filter
END IF
END DO find_j_index_of_fft_filter
grid%mf_fft = grid%msft(ids,j_save)
#endif
CALL pxft ( grid=grid &
,lineno=__LINE__ &
,flag_uv = 0 &
,flag_rurv = 0 &
,flag_wph = 0 &
,flag_ww = 0 &
,flag_t = 1 &
,flag_mu = 0 &
,flag_mut = 0 &
,flag_moist = 0 &
,flag_chem = 0 &
,flag_tracer = 0 &
,flag_scalar = 0 &
,actual_distance_average = .TRUE. &
,pos_def = .FALSE. &
,swap_pole_with_next_j = .FALSE. &
,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
,fft_filter_lat = config_flags%fft_filter_lat &
,dclat = dclat &
,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
grid%ht(i,j) = grid%t_2(i,1,j)
grid%sr(i,j) = grid%sr(i,j) - grid%ht(i,j)
END DO
END DO
#else
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! We stick the topo and map fac in an unused 3d array. The map scale
! factor and computational latitude are passed along for the ride
! (part of the transpose process - we only do 3d arrays) to determine
! "how many" values are used to compute the mean. We want a number
! that is consistent with the original grid resolution.
DO j = jts, MIN(jte,jde-1)
DO k = kts, kte
DO i = its, MIN(ite,ide-1)
grid%t_init(i,k,j) = 1.
END DO
END DO
DO i = its, MIN(ite,ide-1)
grid%t_init(i,1,j) = grid%ht(i,j)
grid%t_init(i,2,j) = grid%msftx(i,j)
grid%t_init(i,3,j) = grid%clat(i,j)
END DO
END DO
# include "XPOSE_POLAR_FILTER_TOPO_z2x.inc"
! Retrieve the 2d arrays for topo, map factors, and the
! computational latitude.
DO j = jpsx, MIN(jpex,jde-1)
DO i = ipsx, MIN(ipex,ide-1)
grid%ht_xxx(i,j) = grid%t_xxx(i,1,j)
grid%mf_xxx(i,j) = grid%t_xxx(i,2,j)
grid%clat_xxx(i,j) = grid%t_xxx(i,3,j)
END DO
END DO
! Get a mean topo field that is consistent with the grid
! distance on each computational latitude loop.
CALL filter_topo ( grid%ht_xxx , grid%clat_xxx , grid%mf_xxx , &
grid%fft_filter_lat , grid%mf_fft , &
.FALSE. , .FALSE. , &
ids, ide, jds, jde, 1 , 1 , &
imsx, imex, jmsx, jmex, 1, 1, &
ipsx, ipex, jpsx, jpex, 1, 1 )
! Stick the filtered topo back into the dummy 3d array to
! transpose it back to "all z on a patch".
DO j = jpsx, MIN(jpex,jde-1)
DO i = ipsx, MIN(ipex,ide-1)
grid%t_xxx(i,1,j) = grid%ht_xxx(i,j)
END DO
END DO
# include "XPOSE_POLAR_FILTER_TOPO_x2z.inc"
! Get the un-transposed topo data.
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
grid%ht(i,j) = grid%t_init(i,1,j)
END DO
END DO
#else
CALL filter_topo ( grid%ht , grid%clat , grid%msftx , &
grid%fft_filter_lat , grid%mf_fft , &
.FALSE. , .FALSE. , &
ids, ide, jds, jde, 1,1, &
ims, ime, jms, jme, 1,1, &
its, ite, jts, jte, 1,1 )
#endif
#endif
ELSE IF ( ( config_flags%map_proj .NE. PROJ_CASSINI ) .AND. ( config_flags%polar ) ) THEN
WRITE ( a_message,* ) 'A global domain (polar = true) requires the Cassini projection'
CALL wrf_error_fatal ( a_message )
END IF
! If we have any input low-res surface pressure, we store it.
IF ( flag_psfc .EQ. 1 ) THEN
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%psfc_gc(i,j) = grid%psfc(i,j)
grid%p_gc(i,1,j) = grid%psfc(i,j)
END DO
END DO
END IF
! If we have the low-resolution surface elevation, stick that in the
! "input" locations of the 3d height. We still have the "hi-res" topo
! stuck in the grid%ht array. The grid%landmask if test is required as some sources
! have ZERO elevation over water (thank you very much).
IF ( flag_soilhgt .EQ. 1) THEN
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
! IF ( grid%landmask(i,j) .GT. 0.5 ) THEN
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid%ght_gc(i,1,j) = grid%toposoil(i,j)
grid%ht_gc(i,j)= grid%toposoil(i,j)
! END IF
END DO
END DO
END IF
! The number of vertical levels in the input data. There is no staggering for
! different variables.
num_metgrid_levels = grid%num_metgrid_levels
! For AFWA UM data, swap incoming extra (theta-based) pressure with the standardly
! named (rho-based) pressure.
IF ( flag_ptheta .EQ. 1 ) THEN
DO j = jts, MIN(jte,jde-1)
DO k = 1 , num_metgrid_levels
DO i = its, MIN(ite,ide-1)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
ptemp = grid%p_gc(i,k,j)
grid%p_gc(i,k,j) = grid%prho_gc(i,k,j)
grid%prho_gc(i,k,j) = ptemp
END DO
END DO
END DO
END IF
! For UM data, the "surface" and the "first hybrid" level for the theta-level data fields are the same.
! Average the surface (k=1) and the second hybrid level (k=num_metgrid_levels-1) to get the first hybrid
! layer. We only do this for the theta-level data: pressure, temperature, specific humidity, and
! geopotential height (i.e. we do not modify u, v, or the rho-based pressure).
IF ( ( flag_ptheta .EQ. 1 ) .OR. ( flag_prho .EQ. 1 ) ) THEN
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid% p_gc(i,num_metgrid_levels,j) = ( grid% p_gc(i,1,j) + grid% p_gc(i,num_metgrid_levels-1,j) ) * 0.5
grid% t_gc(i,num_metgrid_levels,j) = ( grid% t_gc(i,1,j) + grid% t_gc(i,num_metgrid_levels-1,j) ) * 0.5
grid%ght_gc(i,num_metgrid_levels,j) = ( grid%ght_gc(i,1,j) + grid%ght_gc(i,num_metgrid_levels-1,j) ) * 0.5
END DO
END DO
IF ( grid%sh_gc(its,1,jts) .LT. 0 ) THEN
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
grid% sh_gc(i,1,j) = 2. * grid% sh_gc(i,num_metgrid_levels,j) - grid% sh_gc(i,num_metgrid_levels-1,j)
END DO