-
Notifications
You must be signed in to change notification settings - Fork 150
/
Copy pathdyn_comp.F90
1849 lines (1519 loc) · 80.9 KB
/
dyn_comp.F90
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
module dyn_comp
! CAM component interfaces to the MPAS Dynamical Core
use shr_kind_mod, only: r8=>shr_kind_r8
use spmd_utils, only: masterproc, mpicom, npes
use physconst, only: pi, gravit, rair, cpair
use pmgrid, only: plev, plevp
use constituents, only: pcnst, cnst_name, cnst_is_a_water_species, cnst_read_iv
use const_init, only: cnst_init_default
use cam_control_mod, only: initial_run
use cam_initfiles, only: initial_file_get_id, topo_file_get_id
use cam_grid_support, only: cam_grid_id, &
cam_grid_get_latvals, cam_grid_get_lonvals
use cam_map_utils, only: iMap
use inic_analytic, only: analytic_ic_active, dyn_set_inic_col
use dyn_tests_utils, only: vcoord=>vc_height
use cam_history, only: addfld, horiz_only
use string_utils, only: date2yyyymmdd, sec2hms, int2str
use ncdio_atm, only: infld
use pio, only: file_desc_t
use cam_pio_utils, only: clean_iodesc_list
use time_manager, only: get_start_date, get_stop_date, get_run_duration, &
timemgr_get_calendar_cf, get_step_size
use cam_logfile, only: iulog
use cam_abortutils, only: endrun
use mpas_timekeeping, only : MPAS_TimeInterval_type
use cam_mpas_subdriver, only: cam_mpas_global_sum_real
use cam_budget, only: cam_budget_em_snapshot, cam_budget_em_register
use phys_control, only: use_gw_front, use_gw_front_igw
implicit none
private
save
public :: &
dyn_import_t, &
dyn_export_t, &
dyn_readnl, &
dyn_register, &
dyn_init, &
dyn_run, &
dyn_final
! Note that the fields in the import and export states are pointers into the MPAS dycore internal
! data structures. These fields have the order of the vertical and horizontal dimensions swapped
! relative to the CAM convention, as well as having the vertical indices ordered from bottom to top
! of atm. An exception is that the export state contains two fields, pmiddry and pintdry, not managed
! by the MPAS infrastructure. These fields are only used by the physics package and are computed
! in the dp_coupling module.
type dyn_import_t
!
! Number of cells, edges, vertices, and vertical layers in this block
!
integer :: nCells ! Number of cells, including halo cells
integer :: nEdges ! Number of edges, including halo edges
integer :: nVertices ! Number of vertices, including halo vertices
integer :: nVertLevels ! Number of vertical layers
integer :: nCellsSolve ! Number of cells, excluding halo cells
integer :: nEdgesSolve ! Number of edges, excluding halo edges
integer :: nVerticesSolve ! Number of vertices, excluding halo vertices
!
! State that is directly prognosed by the dycore
!
real(r8), dimension(:,:), pointer :: uperp ! Normal velocity at edges [m/s] (nver,nedge)
real(r8), dimension(:,:), pointer :: w ! Vertical velocity [m/s] (nver+1,ncol)
real(r8), dimension(:,:), pointer :: theta_m ! Moist potential temperature [K] (nver,ncol)
real(r8), dimension(:,:), pointer :: rho_zz ! Dry density [kg/m^3]
! divided by d(zeta)/dz (nver,ncol)
real(r8), dimension(:,:,:), pointer :: tracers ! Tracers [kg/kg dry air] (nq,nver,ncol)
!
! Index map between MPAS tracers and CAM constituents
!
integer, dimension(:), pointer :: mpas_from_cam_cnst => null() ! indices into CAM constituent array
!
! Base state variables
!
real(r8), dimension(:,:), pointer :: rho_base ! Base-state dry air density [kg/m^3] (nver,ncol)
real(r8), dimension(:,:), pointer :: theta_base ! Base-state potential temperature [K] (nver,ncol)
!
! Indices of tracers
!
integer :: index_qv ! Index in tracers array of water vapor
! mixing ratio
!
! Invariant -- the vertical coordinate in MPAS-A is a height coordinate
!
real(r8), dimension(:,:), pointer :: zint ! Geometric height [m]
! at layer interfaces (nver+1,ncol)
real(r8), dimension(:,:), pointer :: zz ! Vertical coordinate metric [dimensionless]
! at layer midpoints (nver,ncol)
real(r8), dimension(:), pointer :: fzm ! Interp weight from k layer midpoint to k layer
! interface [dimensionless] (nver)
real(r8), dimension(:), pointer :: fzp ! Interp weight from k-1 layer midpoint to k
! layer interface [dimensionless] (nver)
!
! Invariant -- cell area
!
real(r8), dimension(:), pointer :: areaCell ! cell area (m^2)
!
! Invariant -- needed to compute edge-normal velocities
!
real(r8), dimension(:,:), pointer :: east ! Cartesian components of unit east vector
! at cell centers [dimensionless] (3,ncol)
real(r8), dimension(:,:), pointer :: north ! Cartesian components of unit north vector
! at cell centers [dimensionless] (3,ncol)
real(r8), dimension(:,:), pointer :: normal ! Cartesian components of the vector normal
! to an edge and tangential to the surface
! of the sphere [dimensionless] (3,ncol)
integer, dimension(:,:), pointer :: cellsOnEdge ! Indices of cells separated by an edge (2,nedge)
!
! State that may be directly derived from dycore prognostic state
!
real(r8), dimension(:,:), pointer :: theta ! Potential temperature [K] (nver,ncol)
real(r8), dimension(:,:), pointer :: exner ! Exner function [-] (nver,ncol)
real(r8), dimension(:,:), pointer :: rho ! Dry density [kg/m^3] (nver,ncol)
real(r8), dimension(:,:), pointer :: ux ! Zonal veloc at center [m/s] (nver,ncol)
real(r8), dimension(:,:), pointer :: uy ! Meridional veloc at center [m/s] (nver,ncol)
!
! Tendencies from physics
!
real(r8), dimension(:,:), pointer :: ru_tend ! Normal horizontal momentum tendency
! from physics [kg/m^2/s] (nver,nedge)
real(r8), dimension(:,:), pointer :: rtheta_tend ! Tendency of rho*theta/zz
! from physics [kg K/m^3/s] (nver,ncol)
real(r8), dimension(:,:), pointer :: rho_tend ! Dry air density tendency
! from physics [kg/m^3/s] (nver,ncol)
end type dyn_import_t
type dyn_export_t
!
! Number of cells, edges, vertices, and vertical layers in this block
!
integer :: nCells ! Number of cells, including halo cells
integer :: nEdges ! Number of edges, including halo edges
integer :: nVertices ! Number of vertices, including halo vertices
integer :: nVertLevels ! Number of vertical layers
integer :: nCellsSolve ! Number of cells, excluding halo cells
integer :: nEdgesSolve ! Number of edges, excluding halo edges
integer :: nVerticesSolve ! Number of vertices, excluding halo vertices
!
! State that is directly prognosed by the dycore
!
real(r8), dimension(:,:), pointer :: uperp ! Normal velocity at edges [m/s] (nver,nedge)
real(r8), dimension(:,:), pointer :: w ! Vertical velocity [m/s] (nver+1,ncol)
real(r8), dimension(:,:), pointer :: theta_m ! Moist potential temperature [K] (nver,ncol)
real(r8), dimension(:,:), pointer :: rho_zz ! Dry density [kg/m^3]
! divided by d(zeta)/dz (nver,ncol)
real(r8), dimension(:,:,:), pointer :: tracers ! Tracers [kg/kg dry air] (nq,nver,ncol)
!
! Indices of tracers
!
integer :: index_qv ! Index in tracers array of water vapor
! mixing ratio
!
! Index map between MPAS tracers and CAM constituents
!
integer, dimension(:), pointer :: cam_from_mpas_cnst => null() ! indices into MPAS tracers array
!
! Invariant -- the vertical coordinate in MPAS-A is a height coordinate
!
real(r8), dimension(:,:), pointer :: zint ! Geometric height [m]
! at layer interfaces (nver+1,ncol)
real(r8), dimension(:,:), pointer :: zz ! Vertical coordinate metric [dimensionless]
! at layer midpoints (nver,ncol)
real(r8), dimension(:), pointer :: fzm ! Interp weight from k layer midpoint to k layer
! interface [dimensionless] (nver)
real(r8), dimension(:), pointer :: fzp ! Interp weight from k-1 layer midpoint to k
!
! Invariant -- needed for computing the frontogenesis function
!
real(r8), dimension(:,:), pointer :: defc_a
real(r8), dimension(:,:), pointer :: defc_b
real(r8), dimension(:,:), pointer :: cell_gradient_coef_x
real(r8), dimension(:,:), pointer :: cell_gradient_coef_y
real(r8), dimension(:,:), pointer :: edgesOnCell_sign
real(r8), dimension(:), pointer :: dvEdge
real(r8), dimension(:), pointer :: areaCell ! cell area (m^2)
integer, dimension(:,:), pointer :: edgesOnCell
integer, dimension(:,:), pointer :: cellsOnEdge
integer, dimension(:), pointer :: nEdgesOnCell
real(r8), dimension(:,:), pointer :: utangential ! velocity tangent to cell edge,
! diagnosed by mpas
!
! State that may be directly derived from dycore prognostic state
!
real(r8), dimension(:,:), pointer :: theta ! Potential temperature [K] (nver,ncol)
real(r8), dimension(:,:), pointer :: exner ! Exner function [-] (nver,ncol)
real(r8), dimension(:,:), pointer :: rho ! Dry density [kg/m^3] (nver,ncol)
real(r8), dimension(:,:), pointer :: ux ! Zonal veloc at center [m/s] (nver,ncol)
real(r8), dimension(:,:), pointer :: uy ! Meridional veloc at center [m/s] (nver,ncol)
real(r8), dimension(:,:), pointer :: pmiddry ! Dry hydrostatic pressure [Pa]
! at layer midpoints (nver,ncol)
real(r8), dimension(:,:), pointer :: pintdry ! Dry hydrostatic pressure [Pa]
! at layer interfaces (nver+1,ncol)
real(r8), dimension(:,:), pointer :: vorticity ! Relative vertical vorticity [s^-1]
! (nver,nvtx)
real(r8), dimension(:,:), pointer :: divergence ! Horizontal velocity divergence [s^-1]
! (nver,ncol)
end type dyn_export_t
! Frontogenesis indices
integer, public :: frontgf_idx = -1
integer, public :: frontga_idx = -1
real(r8), parameter :: rad2deg = 180.0_r8 / pi
real(r8), parameter :: deg2rad = pi / 180.0_r8
! The global cell indices are used to seed the RNG which is used to apply
! random perturbations to the initial temperature field. These global indices
! are just those for the local dynamics block.
integer, allocatable :: glob_ind(:)
type (MPAS_TimeInterval_type) :: integrationLength ! set to CAM's dynamics/physics coupling interval
!=========================================================================================
contains
!=========================================================================================
subroutine dyn_readnl(NLFileName)
! Read the dycore-relevant namelists from the input file.
! First must set up basic MPAS infrastructure to allow the MPAS-A dycore
! to save namelist options into MPAS-native datastructures called "pools".
use units, only: getunit
use cam_pio_utils, only: pio_subsystem
use cam_mpas_subdriver, only : domain_ptr, cam_mpas_init_phase1, cam_mpas_init_phase2
use mpas_pool_routines, only : mpas_pool_add_config
! Dummy argument
character(len=*), intent(in) :: NLFileName
! Local variables
integer, dimension(2) :: logUnits ! stdout and stderr for MPAS logging
integer :: yr, mon, day, tod, ndate, nday, nsec
character(len=*), parameter :: subname = 'dyn_comp:dyn_readnl'
!----------------------------------------------------------------------------
logUnits(1) = iulog
logUnits(2) = getunit()
call cam_mpas_init_phase1(mpicom, endrun, logUnits, r8)
! read namelist
call cam_mpas_namelist_read(NLFileName, domain_ptr % configs)
! Set config_start_date, etc. (these will not appear in the dycore namelist)
call get_start_date(yr, mon, day, tod)
ndate = yr*10000 + mon*100 + day
call mpas_pool_add_config(domain_ptr % configs, 'config_start_time', date2yyyymmdd(ndate)//'_'//sec2hms(tod))
call get_stop_date(yr, mon, day, tod)
ndate = yr*10000 + mon*100 + day
call mpas_pool_add_config(domain_ptr % configs, 'config_stop_time', date2yyyymmdd(ndate)//'_'//sec2hms(tod))
call get_run_duration(nday, nsec)
call mpas_pool_add_config(domain_ptr % configs, 'config_run_duration', trim(int2str(nday))//'_'//sec2hms(nsec))
! Although the following namelist options are not expected to be used by CAM-MPAS, the MPAS-A dycore
! references these options, and they therefore must be defined in the configs pool
call mpas_pool_add_config(domain_ptr % configs, 'config_restart_timestamp_name', 'restart_timestamp')
call mpas_pool_add_config(domain_ptr % configs, 'config_IAU_option', 'off')
call mpas_pool_add_config(domain_ptr % configs, 'config_do_DAcycling', .false.)
call mpas_pool_add_config(domain_ptr % configs, 'config_halo_exch_method', 'mpas_halo')
call cam_mpas_init_phase2(pio_subsystem, endrun, timemgr_get_calendar_cf())
end subroutine dyn_readnl
!=========================================================================================
subroutine dyn_register()
! Register fields that are computed by the dycore and passed to the physics via the
! physics buffer.
use physics_buffer, only: pbuf_add_field, dtype_r8
use ppgrid, only: pcols, pver
use phys_control, only: use_gw_front, use_gw_front_igw
!----------------------------------------------------------------------------
! These fields are computed by the dycore and passed to the physics via the
! physics buffer.
if (use_gw_front .or. use_gw_front_igw) then
call pbuf_add_field("FRONTGF", "global", dtype_r8, (/pcols,pver/), frontgf_idx)
call pbuf_add_field("FRONTGA", "global", dtype_r8, (/pcols,pver/), frontga_idx)
end if
end subroutine dyn_register
!=========================================================================================
subroutine dyn_init(dyn_in, dyn_out)
use air_composition, only : thermodynamic_active_species_idx, thermodynamic_active_species_idx_dycore
use air_composition, only : thermodynamic_active_species_num
use air_composition, only : thermodynamic_active_species_liq_idx,thermodynamic_active_species_ice_idx
use air_composition, only : thermodynamic_active_species_liq_idx_dycore,thermodynamic_active_species_ice_idx_dycore
use air_composition, only : thermodynamic_active_species_liq_num, thermodynamic_active_species_ice_num
use cam_mpas_subdriver, only : domain_ptr, cam_mpas_init_phase4
use cam_mpas_subdriver, only : cam_mpas_define_scalars
use mpas_pool_routines, only : mpas_pool_get_subpool, mpas_pool_get_array, mpas_pool_get_dimension, &
mpas_pool_get_config
use mpas_timekeeping, only : MPAS_set_timeInterval
use mpas_derived_types, only : mpas_pool_type
use mpas_constants, only : mpas_constants_compute_derived
use dyn_tests_utils, only : vc_dycore, vc_height, string_vc, vc_str_lgth
use cam_budget, only : thermo_budget_history
! arguments:
type(dyn_import_t), intent(inout) :: dyn_in
type(dyn_export_t), intent(inout) :: dyn_out
! Local variables:
integer :: ierr
type(mpas_pool_type), pointer :: mesh_pool
type(mpas_pool_type), pointer :: state_pool
type(mpas_pool_type), pointer :: diag_pool
type(mpas_pool_type), pointer :: tend_physics_pool
integer, pointer :: nCells
integer, pointer :: nEdges
integer, pointer :: nVertices
integer, pointer :: nVertLevels
integer, pointer :: nCellsSolve
integer, pointer :: nEdgesSolve
integer, pointer :: nVerticesSolve
integer, pointer :: index_qv
integer, pointer :: indexToCellID(:) ! global indices of cell centers of local block
real(r8) :: dtime
real(r8), pointer :: mpas_dt
real(r8) :: dt_ratio
character(len=128) :: errmsg
character(len=*), parameter :: subname = 'dyn_comp::dyn_init'
! variables for initializing energy and axial angular momentum diagnostics
integer, parameter :: num_stages = 6
character (len = 8), dimension(num_stages) :: stage = (/"dBF ","dAP ","dAM ","BD_dparm","BD_DMEA ","BD_phys "/)
character (len = 55),dimension(num_stages) :: stage_txt = (/&
" dynamics state before physics (d_p_coupling) ",&
" dynamics state with T,u,V increment but not q ",&
" dynamics state with full physics increment (incl.q)",&
"dE/dt params+efix in dycore (dparam)(dAP-dBF) ",&
"dE/dt dry mass adjustment in dycore (dAM-dAP)",&
"dE/dt physics total in dycore (phys) (dAM-dBF)" &
/)
integer :: istage, ivars, m
character (len=108) :: str1, str2, str3
character (len=vc_str_lgth) :: vc_str
!-------------------------------------------------------
vc_dycore = vc_height
if (masterproc) then
call string_vc(vc_dycore,vc_str)
write(iulog,*)'vertical coordinate dycore : ',trim(vc_str)
end if
!----------------------------------------------------------------------------
if (initial_run) then
call cam_mpas_define_scalars(domain_ptr % blocklist, dyn_in % mpas_from_cam_cnst, &
dyn_out % cam_from_mpas_cnst, ierr)
if (ierr /= 0) then
call endrun(subname//': Set-up of constituents for MPAS-A dycore failed.')
end if
end if
call mpas_pool_get_subpool(domain_ptr % blocklist % structs, 'mesh', mesh_pool)
call mpas_pool_get_subpool(domain_ptr % blocklist % structs, 'state', state_pool)
call mpas_pool_get_subpool(domain_ptr % blocklist % structs, 'diag', diag_pool)
call mpas_pool_get_subpool(domain_ptr % blocklist % structs, 'tend_physics', tend_physics_pool)
! Let dynamics import state point to memory managed by MPAS-Atmosphere
call mpas_pool_get_dimension(mesh_pool, 'nCells', nCells)
dyn_in % nCells = nCells
call mpas_pool_get_dimension(mesh_pool, 'nEdges', nEdges)
dyn_in % nEdges = nEdges
call mpas_pool_get_dimension(mesh_pool, 'nVertices', nVertices)
dyn_in % nVertices = nVertices
call mpas_pool_get_dimension(mesh_pool, 'nVertLevels', nVertLevels)
dyn_in % nVertLevels = nVertLevels
call mpas_pool_get_dimension(mesh_pool, 'nCellsSolve', nCellsSolve)
dyn_in % nCellsSolve = nCellsSolve
call mpas_pool_get_dimension(mesh_pool, 'nEdgesSolve', nEdgesSolve)
dyn_in % nEdgesSolve = nEdgesSolve
call mpas_pool_get_dimension(mesh_pool, 'nVerticesSolve', nVerticesSolve)
dyn_in % nVerticesSolve = nVerticesSolve
! In MPAS timeLevel=1 is the current state. So the fields input to the dycore should
! be in timeLevel=1.
call mpas_pool_get_array(state_pool, 'u', dyn_in % uperp, timeLevel=1)
call mpas_pool_get_array(state_pool, 'w', dyn_in % w, timeLevel=1)
call mpas_pool_get_array(state_pool, 'theta_m', dyn_in % theta_m, timeLevel=1)
call mpas_pool_get_array(state_pool, 'rho_zz', dyn_in % rho_zz, timeLevel=1)
call mpas_pool_get_array(state_pool, 'scalars', dyn_in % tracers, timeLevel=1)
call mpas_pool_get_array(diag_pool, 'rho_base', dyn_in % rho_base)
call mpas_pool_get_array(diag_pool, 'theta_base', dyn_in % theta_base)
call mpas_pool_get_dimension(state_pool, 'index_qv', index_qv)
dyn_in % index_qv = index_qv
call mpas_pool_get_array(mesh_pool, 'zgrid', dyn_in % zint)
call mpas_pool_get_array(mesh_pool, 'zz', dyn_in % zz)
call mpas_pool_get_array(mesh_pool, 'fzm', dyn_in % fzm)
call mpas_pool_get_array(mesh_pool, 'fzp', dyn_in % fzp)
call mpas_pool_get_array(mesh_pool, 'areaCell', dyn_in % areaCell)
call mpas_pool_get_array(mesh_pool, 'east', dyn_in % east)
call mpas_pool_get_array(mesh_pool, 'north', dyn_in % north)
call mpas_pool_get_array(mesh_pool, 'edgeNormalVectors', dyn_in % normal)
call mpas_pool_get_array(mesh_pool, 'cellsOnEdge', dyn_in % cellsOnEdge)
call mpas_pool_get_array(diag_pool, 'theta', dyn_in % theta)
call mpas_pool_get_array(diag_pool, 'exner', dyn_in % exner)
call mpas_pool_get_array(diag_pool, 'rho', dyn_in % rho)
call mpas_pool_get_array(diag_pool, 'uReconstructZonal', dyn_in % ux)
call mpas_pool_get_array(diag_pool, 'uReconstructMeridional', dyn_in % uy)
! Let dynamics export state point to memory managed by MPAS-Atmosphere
! Exception: pmiddry and pintdry are not managed by the MPAS infrastructure
dyn_out % nCells = dyn_in % nCells
dyn_out % nEdges = dyn_in % nEdges
dyn_out % nVertices = dyn_in % nVertices
dyn_out % nVertLevels = dyn_in % nVertLevels
dyn_out % nCellsSolve = dyn_in % nCellsSolve
dyn_out % nEdgesSolve = dyn_in % nEdgesSolve
dyn_out % nVerticesSolve = dyn_in % nVerticesSolve
dyn_out % index_qv = dyn_in % index_qv
! MPAS swaps pointers internally so that after a dycore timestep, the updated state is
! in timeLevel=1. Thus we want dyn_out to also point to timeLevel=1. Can just copy
! the pointers from dyn_in.
dyn_out % uperp => dyn_in % uperp
dyn_out % w => dyn_in % w
dyn_out % theta_m => dyn_in % theta_m
dyn_out % rho_zz => dyn_in % rho_zz
dyn_out % tracers => dyn_in % tracers
! These components don't have a time level index.
dyn_out % zint => dyn_in % zint
dyn_out % zz => dyn_in % zz
dyn_out % fzm => dyn_in % fzm
dyn_out % fzp => dyn_in % fzp
dyn_out % theta => dyn_in % theta
dyn_out % exner => dyn_in % exner
dyn_out % rho => dyn_in % rho
dyn_out % ux => dyn_in % ux
dyn_out % uy => dyn_in % uy
! for frontogenesis calc
if (use_gw_front .or. use_gw_front_igw) then
dyn_out % areaCell => dyn_in % areaCell
dyn_out % cellsOnEdge => dyn_in % cellsOnEdge
call mpas_pool_get_array(mesh_pool, 'defc_a', dyn_out % defc_a)
call mpas_pool_get_array(mesh_pool, 'defc_b', dyn_out % defc_b)
call mpas_pool_get_array(mesh_pool, 'cell_gradient_coef_x', dyn_out % cell_gradient_coef_x)
call mpas_pool_get_array(mesh_pool, 'cell_gradient_coef_y', dyn_out % cell_gradient_coef_y)
call mpas_pool_get_array(mesh_pool, 'edgesOnCell_sign', dyn_out % edgesOnCell_sign)
call mpas_pool_get_array(mesh_pool, 'dvEdge', dyn_out % dvEdge)
call mpas_pool_get_array(mesh_pool, 'edgesOnCell', dyn_out % edgesOnCell)
call mpas_pool_get_array(mesh_pool, 'nEdgesOnCell', dyn_out % nEdgesOnCell)
call mpas_pool_get_array(diag_pool, 'v', dyn_out % utangential)
endif
! cam-required hydrostatic pressures
allocate(dyn_out % pmiddry(nVertLevels, nCells), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate dyn_out%pmiddry array')
allocate(dyn_out % pintdry(nVertLevels+1, nCells), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate dyn_out%pintdry array')
call mpas_pool_get_array(diag_pool, 'vorticity', dyn_out % vorticity)
call mpas_pool_get_array(diag_pool, 'divergence', dyn_out % divergence)
call mpas_pool_get_array(mesh_pool, 'indexToCellID', indexToCellID)
allocate(glob_ind(nCellsSolve), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate glob_ind array')
glob_ind = indexToCellID(1:nCellsSolve)
call mpas_constants_compute_derived()
if (initial_run) then
call read_inidat(dyn_in)
call clean_iodesc_list()
end if
call cam_mpas_init_phase4(endrun)
!
! Set pointers to tendency fields that are not allocated until the call to cam_mpas_init_phase4
!
call mpas_pool_get_array(tend_physics_pool, 'tend_ru_physics', dyn_in % ru_tend)
call mpas_pool_get_array(tend_physics_pool, 'tend_rtheta_physics', dyn_in % rtheta_tend)
call mpas_pool_get_array(tend_physics_pool, 'tend_rho_physics', dyn_in % rho_tend)
! Check that CAM's timestep, i.e., the dynamics/physics coupling interval, is an integer multiple
! of the MPAS timestep.
! Get CAM time step
dtime = get_step_size()
! Get MPAS-A dycore time step
call mpas_pool_get_config(domain_ptr % configs, 'config_dt', mpas_dt)
! Calculate time step ratio
dt_ratio = dtime / mpas_dt
! Stop if the dycore time step does not evenly divide the CAM time step
if (ceiling(dt_ratio) /= floor(dt_ratio)) then
write(errmsg, '(a,f9.3,a,f9.3,a)') 'The ratio of the CAM timestep, ', dtime, &
' to the MPAS-A dycore timestep, ', mpas_dt, ' is not an integer'
call endrun(subname//': '//trim(errmsg))
end if
! dtime has no fractional part, but use nint to deal with any roundoff errors.
! Set the interval over which the dycore should integrate during each call to dyn_run.
call MPAS_set_timeInterval(integrationLength, S=nint(dtime), S_n=0, S_d=1)
!
! initialize history for MPAS energy budgets
if (thermo_budget_history) then
! Define energy/mass snapshots using stage structure
do istage = 1, num_stages
call cam_budget_em_snapshot(TRIM(ADJUSTL(stage(istage))), 'dyn', longname=TRIM(ADJUSTL(stage_txt(istage))))
end do
!
! initialize MPAS energy budgets
! add budgets that are derived from stages
!
call cam_budget_em_register('dEdt_param_efix_in_dyn','dAP','dBF',pkgtype='dyn',optype='dif', &
longname="dE/dt parameterizations+efix in dycore (dparam)(dAP-dBF)")
call cam_budget_em_register('dEdt_dme_adjust_in_dyn','dAM','dAP',pkgtype='dyn',optype='dif', &
longname="dE/dt dry mass adjustment in dycore (dAM-dAP)")
call cam_budget_em_register('dEdt_phys_total_in_dyn','dAM','dBF',pkgtype='dyn',optype='dif', &
longname="dE/dt physics total in dycore (phys) (dAM-dBF)")
end if
!
! initialize CAM thermodynamic infrastructure
!
do m=1,thermodynamic_active_species_num
thermodynamic_active_species_idx_dycore(m) = dyn_out % cam_from_mpas_cnst(thermodynamic_active_species_idx(m))
if (masterproc) then
write(iulog,'(a,2I4)') subname//": m,thermodynamic_active_species_idx_dycore: ", &
m,thermodynamic_active_species_idx_dycore(m)
end if
end do
do m=1,thermodynamic_active_species_liq_num
thermodynamic_active_species_liq_idx_dycore(m) = dyn_out % cam_from_mpas_cnst(thermodynamic_active_species_liq_idx(m))
if (masterproc) then
write(iulog,'(a,2I4)') subname//": m,thermodynamic_active_species_idx_liq_dycore: ", &
m,thermodynamic_active_species_liq_idx_dycore(m)
end if
end do
do m=1,thermodynamic_active_species_ice_num
thermodynamic_active_species_ice_idx_dycore(m) = dyn_out % cam_from_mpas_cnst(thermodynamic_active_species_ice_idx(m))
if (masterproc) then
write(iulog,'(a,2I4)') subname//": m,thermodynamic_active_species_idx_ice_dycore: ", &
m,thermodynamic_active_species_ice_idx_dycore(m)
end if
end do
end subroutine dyn_init
!=========================================================================================
subroutine dyn_run(dyn_in, dyn_out)
use cam_mpas_subdriver, only : cam_mpas_run
use cam_mpas_subdriver, only : domain_ptr
use mpas_pool_routines, only : mpas_pool_get_subpool, mpas_pool_get_array
use mpas_derived_types, only : mpas_pool_type
! Advances the dynamics state provided in dyn_in by one physics
! timestep to produce dynamics state held in dyn_out.
type (dyn_import_t), intent(inout) :: dyn_in
type (dyn_export_t), intent(inout) :: dyn_out
! Local variables
type(mpas_pool_type), pointer :: state_pool
character(len=*), parameter :: subname = 'dyn_comp:dyn_run'
real(r8) :: dtime
!----------------------------------------------------------------------------
! Call the MPAS-A dycore
call cam_mpas_run(integrationLength)
! Update the dyn_in/dyn_out pointers to the current state of the prognostic fields.
call mpas_pool_get_subpool(domain_ptr % blocklist % structs, 'state', state_pool)
call mpas_pool_get_array(state_pool, 'u', dyn_in % uperp, timeLevel=1)
call mpas_pool_get_array(state_pool, 'w', dyn_in % w, timeLevel=1)
call mpas_pool_get_array(state_pool, 'theta_m', dyn_in % theta_m, timeLevel=1)
call mpas_pool_get_array(state_pool, 'rho_zz', dyn_in % rho_zz, timeLevel=1)
call mpas_pool_get_array(state_pool, 'scalars', dyn_in % tracers, timeLevel=1)
dyn_out % uperp => dyn_in % uperp
dyn_out % w => dyn_in % w
dyn_out % theta_m => dyn_in % theta_m
dyn_out % rho_zz => dyn_in % rho_zz
dyn_out % tracers => dyn_in % tracers
end subroutine dyn_run
subroutine dyn_final(dyn_in, dyn_out)
use cam_mpas_subdriver, only : cam_mpas_finalize
! Deallocates the dynamics import and export states, and finalizes
! the MPAS dycore.
type (dyn_import_t), intent(inout) :: dyn_in
type (dyn_export_t), intent(inout) :: dyn_out
character(len=*), parameter :: subname = 'dyn_comp:dyn_final'
!----------------------------------------------------------------------------
!
! Prevent any further access to MPAS-Atmosphere memory
!
dyn_in % nCells = 0
dyn_in % nEdges = 0
dyn_in % nVertices = 0
dyn_in % nVertLevels = 0
dyn_in % nCellsSolve = 0
dyn_in % nEdgesSolve = 0
dyn_in % nVerticesSolve = 0
nullify(dyn_in % uperp)
nullify(dyn_in % w)
nullify(dyn_in % theta_m)
nullify(dyn_in % rho_zz)
nullify(dyn_in % tracers)
deallocate(dyn_in % mpas_from_cam_cnst)
nullify(dyn_in % rho_base)
nullify(dyn_in % theta_base)
dyn_in % index_qv = 0
nullify(dyn_in % zint)
nullify(dyn_in % zz)
nullify(dyn_in % fzm)
nullify(dyn_in % fzp)
nullify(dyn_in % east)
nullify(dyn_in % north)
nullify(dyn_in % normal)
nullify(dyn_in % cellsOnEdge)
nullify(dyn_in % theta)
nullify(dyn_in % exner)
nullify(dyn_in % rho)
nullify(dyn_in % ux)
nullify(dyn_in % uy)
nullify(dyn_in % ru_tend)
nullify(dyn_in % rtheta_tend)
nullify(dyn_in % rho_tend)
!
! Prevent any further access to MPAS-Atmosphere memory
!
dyn_out % nCells = 0
dyn_out % nEdges = 0
dyn_out % nVertices = 0
dyn_out % nVertLevels = 0
dyn_out % nCellsSolve = 0
dyn_out % nEdgesSolve = 0
dyn_out % nVerticesSolve = 0
nullify(dyn_out % uperp)
nullify(dyn_out % w)
nullify(dyn_out % theta_m)
nullify(dyn_out % rho_zz)
nullify(dyn_out % tracers)
deallocate(dyn_out % cam_from_mpas_cnst)
dyn_out % index_qv = 0
nullify(dyn_out % zint)
nullify(dyn_out % zz)
nullify(dyn_out % fzm)
nullify(dyn_out % fzp)
nullify(dyn_out % theta)
nullify(dyn_out % exner)
nullify(dyn_out % rho)
nullify(dyn_out % ux)
nullify(dyn_out % uy)
deallocate(dyn_out % pmiddry)
deallocate(dyn_out % pintdry)
nullify(dyn_out % vorticity)
nullify(dyn_out % divergence)
call cam_mpas_finalize()
end subroutine dyn_final
!=========================================================================================
! Private routines.
!=========================================================================================
subroutine read_inidat(dyn_in)
! Set initial conditions. Either from analytic expressions or read from file.
use cam_mpas_subdriver, only : domain_ptr, cam_mpas_update_halo, cam_mpas_cell_to_edge_winds
use cam_initfiles, only : scale_dry_air_mass
use mpas_pool_routines, only : mpas_pool_get_subpool, mpas_pool_get_array
use mpas_derived_types, only : mpas_pool_type
use mpas_vector_reconstruction, only : mpas_reconstruct
use mpas_constants, only : Rv_over_Rd => rvord
use mpas_constants, only : rgas
use mpas_constants, only : p0
use mpas_constants, only : gravity
use string_utils, only : int2str
use shr_kind_mod, only : shr_kind_cx
! arguments
type(dyn_import_t), target, intent(inout) :: dyn_in
! Local variables
integer :: nCellsSolve, nEdgesSolve
integer :: i, k, kk, m
type(file_desc_t), pointer :: fh_ini
type(file_desc_t), pointer :: fh_topo
real(r8), allocatable :: latvals(:)
real(r8), allocatable :: lonvals(:)
real(r8), pointer :: latvals_deg(:)
real(r8), pointer :: lonvals_deg(:)
real(r8), pointer :: uperp(:,:) ! Normal velocity at edges [m/s] (nver,nedge)
real(r8), pointer :: w(:,:) ! Vertical velocity [m/s] (nver+1,ncol)
real(r8), pointer :: theta_m(:,:) ! Moist potential temperature [K] (nver,ncol)
real(r8), pointer :: rho_zz(:,:) ! Dry density [kg/m^3]
! divided by d(zeta)/dz (nver,ncol)
real(r8), pointer :: tracers(:,:,:) ! Tracers [kg/kg dry air] (nq,nver,ncol)
real(r8), pointer :: zint(:,:) ! Geometric height [m]
! at layer interfaces (nver+1,ncol)
real(r8), pointer :: zz(:,:) ! Vertical coordinate metric [dimensionless]
! at layer midpoints (nver,ncol)
real(r8), pointer :: theta(:,:) ! Potential temperature [K] (nver,ncol)
real(r8), pointer :: rho(:,:) ! Dry density [kg/m^3] (nver,ncol)
real(r8), pointer :: ux(:,:) ! Zonal veloc at center [m/s] (nver,ncol)
real(r8), pointer :: uy(:,:) ! Meridional veloc at center [m/s] (nver,ncol)
real(r8), pointer :: theta_base(:,:)
real(r8), pointer :: rho_base(:,:)
integer :: ixqv
integer, dimension(:), pointer :: mpas_from_cam_cnst
integer, allocatable :: m_ind(:)
real(r8), allocatable :: &
cam2d(:), cam3d(:,:), cam4d(:,:,:), zi(:,:) ! temp arrays using CAM data order
real(r8), allocatable :: zsurf(:)
! temp arrays using MPAS data order
real(r8), allocatable :: t(:,:) ! temperature
real(r8), allocatable :: pintdry(:,:) ! dry interface pressures
real(r8), allocatable :: pmiddry(:,:) ! dry midpoint pressures
real(r8), allocatable :: pmid(:,:) ! midpoint pressures
real(r8), allocatable :: mpas3d(:,:,:)
real(r8), allocatable :: qv(:), tm(:)
logical :: readvar
character(len=shr_kind_cx) :: str
type(mpas_pool_type), pointer :: mesh_pool
type(mpas_pool_type), pointer :: diag_pool
real(r8), pointer :: uReconstructX(:,:)
real(r8), pointer :: uReconstructY(:,:)
real(r8), pointer :: uReconstructZ(:,:)
integer :: mpas_idx, cam_idx, ierr
character(len=32) :: trac_name
character(len=*), parameter :: subname = 'dyn_comp:read_inidat'
!--------------------------------------------------------------------------------------
fh_ini => initial_file_get_id()
fh_topo => topo_file_get_id()
nCellsSolve = dyn_in % nCellsSolve
nEdgesSolve = dyn_in % nEdgesSolve
ixqv = dyn_in % index_qv
mpas_from_cam_cnst => dyn_in % mpas_from_cam_cnst
uperp => dyn_in % uperp
w => dyn_in % w
theta_m => dyn_in % theta_m
rho_zz => dyn_in % rho_zz
tracers => dyn_in % tracers
zint => dyn_in % zint
zz => dyn_in % zz
theta => dyn_in % theta
rho => dyn_in % rho
ux => dyn_in % ux
uy => dyn_in % uy
rho_base => dyn_in % rho_base
theta_base => dyn_in % theta_base
! Check that number of advected tracers is consistent with MPAS.
if (pcnst /= size(tracers, 1)) then
write(iulog,*) subname//': number of tracers, pcnst:', size(tracers,1), pcnst
call endrun(subname//': number of tracers /= pcnst')
end if
! lat/lon needed in radians
latvals_deg => cam_grid_get_latvals(cam_grid_id('mpas_cell'))
lonvals_deg => cam_grid_get_lonvals(cam_grid_id('mpas_cell'))
allocate(latvals(nCellsSolve), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate latvals array')
allocate(lonvals(nCellsSolve), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate lonvals array')
latvals(:) = latvals_deg(:)*deg2rad
lonvals(:) = lonvals_deg(:)*deg2rad
! Set ICs. Either from analytic expressions or read from file.
allocate( &
! temporary arrays using CAM indexing
cam2d(nCellsSolve), &
cam3d(nCellsSolve,plev), &
cam4d(nCellsSolve,plev,pcnst), &
zi(nCellsSolve,plevp), &
! temporary arrays using MPAS indexing
t(plev,nCellsSolve), &
pintdry(plevp,nCellsSolve), &
pmiddry(plev,nCellsSolve), &
pmid(plev,nCellsSolve), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate tmp arrays using CAM and MPAS indexing')
do k = 1, plevp
kk = plevp - k + 1
zi(:,kk) = zint(k,:nCellsSolve)
end do
! If using a topo file check that PHIS is consistent with the surface z coordinate.
if (associated(fh_topo)) then
allocate(zsurf(nCellsSolve), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate zsurf array')
call get_zsurf_from_topo(fh_topo, zsurf)
do i = 1, nCellsSolve
if (abs(zi(i,plevp) - zsurf(i)) > 0.001_r8) then
write(str,*) 'zi= ', zi(i,plevp), ' zsurf= ', zsurf(i),' i= ',i
write(iulog,*) subname//': ERROR: '//TRIM(str)
call endrun(subname//': ERROR: PHIS not consistent with surface z coordinate; '//TRIM(str))
end if
end do
deallocate(zsurf)
else
do i = 1, nCellsSolve
if (abs(zi(i,plevp)) > 1.0E-12_r8) then
write(str,*) 'zi= ', zi(i,plevp), ' but PHIS should be zero'
write(iulog,*) subname//': ERROR: '//TRIM(str)
call endrun(subname//': ERROR: PHIS not consistent with surface z coordinate; '//TRIM(str))
end if
end do
end if
if (analytic_ic_active()) then
w(:,1:nCellsSolve) = 0.0_r8
! U, V cell center velocity components
call dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint=zi, U=cam3d)
do k = 1, plev
kk = plev - k + 1
do i = 1, nCellsSolve
ux(kk,i) = cam3d(i,k)
end do
end do
call dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint=zi, V=cam3d)
do k = 1, plev
kk = plev - k + 1
do i = 1, nCellsSolve
uy(kk,i) = cam3d(i,k)
end do
end do
! Compute uperp by projecting ux and uy from cell centers to edges
call cam_mpas_update_halo('uReconstructZonal', endrun) ! ux => uReconstructZonal
call cam_mpas_update_halo('uReconstructMeridional', endrun) ! uy => uReconstructMeridional
call cam_mpas_cell_to_edge_winds(dyn_in % nEdges, ux, uy, dyn_in % east, dyn_in % north, &
dyn_in % normal, dyn_in % cellsOnEdge, uperp)
call cam_mpas_update_halo('u', endrun) ! u is the name of uperp in the MPAS state pool
! Constituents
allocate(m_ind(pcnst), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate m_ind array')
do m = 1, pcnst
m_ind(m) = m
end do
call dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint=zi, m_cnst=m_ind, Q=cam4d)
do m = 1, pcnst ! index into MPAS tracers array
do k = 1, plev
kk = plev - k + 1
do i = 1, nCellsSolve
tracers(m,kk,i) = cam4d(i,k,mpas_from_cam_cnst(m))
end do
end do
end do
deallocate(m_ind)
! Temperature
call dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint=zi, T=cam3d)
do k = 1, plev
kk = plev - k + 1
do i = 1, nCellsSolve
t(kk,i) = cam3d(i,k)
end do
end do
! Pressures are needed to convert temperature to potential temperature.
call dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint=zi, PS=cam2d)
do i = 1, nCellsSolve
pintdry(1,i) = cam2d(i)
end do
allocate(qv(plev), tm(plev), stat=ierr)
if( ierr /= 0 ) call endrun(subname//': failed to allocate qv and tm arrays')
do i = 1, nCellsSolve
do k = 1, plev
! convert specific humidity to mixing ratio relative to dry air
tracers(1,k,i) = tracers(1,k,i)/(1.0_r8 - tracers(1,k,i))
qv(k) = tracers(1,k,i)
! convert temperature to tm (we are using dry air density and Rd in state eqn)
tm(k) = t(k,i)*(1.0_r8+(Rv_over_Rd)*qv(k))
end do
! integrate up from surface to first mid level. This is full mid-level pressure (i.e. accounts for vapor).
! we are assuming that pintdry(1,i) is the full surface pressure here.
pmid(1,i) = pintdry(1,i)/(1.0_r8+0.5_r8*(zint(2,i)-zint(1,i))*(1.0_r8+qv(1))*gravity/(rgas*tm(1)))
! integrate up the column