-
Notifications
You must be signed in to change notification settings - Fork 0
/
dyn_comp.F90
2127 lines (1811 loc) · 79.8 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 interfaces to the SE Dynamical Core
use shr_kind_mod, only: r8=>shr_kind_r8, shr_kind_cl
use physconst, only: pi
use spmd_utils, only: iam, masterproc
use constituents, only: pcnst, cnst_get_ind, cnst_name, cnst_longname, &
cnst_read_iv, qmin, cnst_type
use cam_control_mod, only: initial_run
use cam_initfiles, only: initial_file_get_id, topo_file_get_id, pertlim
use phys_control, only: use_gw_front, use_gw_front_igw
use dyn_grid, only: timelevel, hvcoord, edgebuf
use cam_grid_support, only: cam_grid_id, cam_grid_get_gcid, &
cam_grid_dimensions, cam_grid_get_dim_names, &
cam_grid_get_latvals, cam_grid_get_lonvals, &
max_hcoordname_len
use cam_map_utils, only: iMap
use inic_analytic, only: analytic_ic_active, analytic_ic_set_ic
use dyn_tests_utils, only: vcoord=>vc_dry_pressure
use cam_history, only: outfld, hist_fld_active, fieldname_len
use cam_history_support, only: max_fieldname_len
use time_manager, only: get_step_size
use ncdio_atm, only: infld
use pio, only: file_desc_t, pio_seterrorhandling, PIO_BCAST_ERROR, &
pio_inq_dimid, pio_inq_dimlen, PIO_NOERR
use infnan, only: isnan
use cam_logfile, only: iulog
use cam_abortutils, only: endrun
use shr_sys_mod, only: shr_sys_flush
use parallel_mod, only: par
use hybrid_mod, only: hybrid_t
use dimensions_mod, only: nelemd, nlev, np, npsq, ntrac, nc, fv_nphys, &
qsize
use element_mod, only: element_t, elem_state_t
use fvm_control_volume_mod, only: fvm_struct, n0_fvm
use time_mod, only: nsplit
use edge_mod, only: initEdgeBuffer, edgeVpack, edgeVunpack, FreeEdgeBuffer
use edgetype_mod, only: EdgeBuffer_t
use bndry_mod, only: bndry_exchange
implicit none
private
save
public :: &
dyn_import_t, &
dyn_export_t, &
dyn_readnl, &
dyn_register, &
dyn_init, &
dyn_run, &
dyn_final
type dyn_import_t
type (element_t), pointer :: elem(:) => null()
type (fvm_struct), pointer :: fvm(:) => null()
end type dyn_import_t
type dyn_export_t
type (element_t), pointer :: elem(:) => null()
type (fvm_struct), pointer :: fvm(:) => null()
end type dyn_export_t
! Namelist
logical, public, protected :: write_restart_unstruct
! Frontogenesis indices
integer, public :: frontgf_idx = -1
integer, public :: frontga_idx = -1
interface read_dyn_var
module procedure read_dyn_field_2d
module procedure read_dyn_field_3d
end interface read_dyn_var
real(r8), parameter :: rad2deg = 180.0_r8 / pi
real(r8), parameter :: deg2rad = pi / 180.0_r8
!===============================================================================
contains
!===============================================================================
subroutine dyn_readnl(NLFileName)
use namelist_utils, only: find_group_name
use namelist_mod, only: homme_set_defaults, homme_postprocess_namelist
use units, only: getunit, freeunit
use spmd_utils, only: masterproc, masterprocid, mpicom, npes
use spmd_utils, only: mpi_real8, mpi_integer, mpi_character, mpi_logical
use dyn_grid, only: se_write_grid_file, se_grid_filename, se_write_gll_corners
use dp_mapping, only: nphys_pts
use native_mapping, only: native_mapping_readnl
use control_mod, only: TRACERTRANSPORT_SE_GLL, tracer_transport_type
use control_mod, only: TRACERTRANSPORT_CONSISTENT_SE_FVM
use control_mod, only: hypervis_subcycle
use control_mod, only: hypervis_subcycle_q, statefreq, runtype
use control_mod, only: nu, nu_div, nu_p, nu_q, nu_top, qsplit, rsplit
use control_mod, only: vert_remap_q_alg, tstep_type, rk_stage_user
use control_mod, only: ftype, limiter_option, partmethod
use control_mod, only: topology, tasknum
use control_mod, only: remap_type
use control_mod, only: fine_ne, hypervis_power, hypervis_scaling
use control_mod, only: max_hypervis_courant, statediag_numtrac,refined_mesh
use control_mod, only: se_met_nudge_u, se_met_nudge_p, se_met_nudge_t, se_met_tevolve
use dimensions_mod, only: qsize_d, ne, npart
use dimensions_mod, only: qsize_condensate_loading, lcp_moist
use dimensions_mod, only: hypervis_on_plevs
use params_mod, only: SFCURVE
use parallel_mod, only: initmpi
use thread_mod, only: initomp, max_num_threads
use thread_mod, only: horz_num_threads, vert_num_threads, tracer_num_threads
! Dummy argument
character(len=*), intent(in) :: NLFileName
! Local variables
integer :: unitn, ierr
real(r8) :: uniform_res_hypervis_scaling,nu_fac
! SE Namelist variables
integer :: se_qsize_condensate_loading
integer :: se_fine_ne
integer :: se_ftype
integer :: se_statediag_numtrac
integer :: se_fv_nphys
real(r8) :: se_hypervis_power
real(r8) :: se_hypervis_scaling
integer :: se_hypervis_subcycle
integer :: se_hypervis_subcycle_q
integer :: se_limiter_option
real(r8) :: se_max_hypervis_courant
character(len=SHR_KIND_CL) :: se_mesh_file
integer :: se_ne
integer :: se_npes
integer :: se_nsplit
real(r8) :: se_nu
real(r8) :: se_nu_div
real(r8) :: se_nu_p
real(r8) :: se_nu_top
integer :: se_qsplit
logical :: se_refined_mesh
integer :: se_rsplit
integer :: se_statefreq
integer :: se_tstep_type
integer :: se_vert_remap_q_alg
integer :: se_horz_num_threads
integer :: se_vert_num_threads
integer :: se_tracer_num_threads
logical :: se_hypervis_on_plevs
logical :: se_write_restart_unstruct
namelist /dyn_se_inparm/ &
se_qsize_condensate_loading, &
se_fine_ne, & ! For refined meshes
se_ftype, & ! forcing type
se_statediag_numtrac, &
se_fv_nphys, &
se_hypervis_power, &
se_hypervis_scaling, &
se_hypervis_subcycle, &
se_hypervis_subcycle_q, &
se_limiter_option, &
se_max_hypervis_courant, &
se_mesh_file, & ! Refined mesh definition file
se_ne, &
se_npes, &
se_nsplit, & ! # of dyn steps per physics timestep
se_nu, &
se_nu_div, &
se_nu_p, &
se_nu_top, &
se_qsplit, &
se_refined_mesh, &
se_rsplit, &
se_statefreq, & ! number of steps per printstate call
se_tstep_type, &
se_vert_remap_q_alg, &
se_met_nudge_u, &
se_met_nudge_p, &
se_met_nudge_t, &
se_met_tevolve, &
se_write_grid_file, &
se_grid_filename, &
se_write_gll_corners, &
se_horz_num_threads, &
se_vert_num_threads, &
se_tracer_num_threads, &
se_hypervis_on_plevs, &
se_write_restart_unstruct
!--------------------------------------------------------------------------
! defaults for variables not set by build-namelist
se_fine_ne = -1
se_hypervis_power = 0
se_hypervis_scaling = 0
se_max_hypervis_courant = 1.0e99_r8
se_mesh_file = ''
se_npes = npes
se_write_restart_unstruct = .false.
! Read the namelist (dyn_se_inparm)
call MPI_barrier(mpicom, ierr)
if (masterproc) then
write(iulog, *) "dyn_readnl: reading dyn_se_inparm namelist..."
unitn = getunit()
open( unitn, file=trim(NLFileName), status='old' )
call find_group_name(unitn, 'dyn_se_inparm', status=ierr)
if (ierr == 0) then
read(unitn, dyn_se_inparm, iostat=ierr)
if (ierr /= 0) then
call endrun('dyn_readnl: ERROR reading dyn_se_inparm namelist')
end if
end if
close(unitn)
call freeunit(unitn)
end if
! Broadcast namelist values to all PEs
call MPI_bcast(se_qsize_condensate_loading, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_fine_ne, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_ftype, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_statediag_numtrac, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_hypervis_power, 1, mpi_real8, masterprocid, mpicom, ierr)
call MPI_bcast(se_hypervis_scaling, 1, mpi_real8, masterprocid, mpicom, ierr)
call MPI_bcast(se_hypervis_subcycle, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_hypervis_subcycle_q, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_limiter_option, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_max_hypervis_courant, 1, mpi_real8, masterprocid, mpicom, ierr)
call MPI_bcast(se_mesh_file, SHR_KIND_CL, mpi_character, masterprocid, mpicom, ierr)
call MPI_bcast(se_ne, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_npes, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_nsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_nu, 1, mpi_real8, masterprocid, mpicom, ierr)
call MPI_bcast(se_nu_div, 1, mpi_real8, masterprocid, mpicom, ierr)
call MPI_bcast(se_nu_p, 1, mpi_real8, masterprocid, mpicom, ierr)
call MPI_bcast(se_nu_top, 1, mpi_real8, masterprocid, mpicom, ierr)
call MPI_bcast(se_qsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_refined_mesh, 1, mpi_logical, masterprocid, mpicom, ierr)
call MPI_bcast(se_rsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_statefreq, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_tstep_type, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_vert_remap_q_alg, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_met_nudge_u, 1, MPI_real8, masterprocid, mpicom,ierr)
call MPI_bcast(se_met_nudge_p, 1, MPI_real8, masterprocid, mpicom,ierr)
call MPI_bcast(se_met_nudge_t, 1, MPI_real8, masterprocid, mpicom,ierr)
call MPI_bcast(se_met_tevolve, 1, MPI_integer, masterprocid, mpicom,ierr)
call MPI_bcast(se_fv_nphys, 1, mpi_integer, masterprocid, mpicom, ierr)
call MPI_bcast(se_write_grid_file, 16, mpi_character, masterprocid, mpicom, ierr)
call MPI_bcast(se_grid_filename, shr_kind_cl, mpi_character, masterprocid, mpicom, ierr)
call MPI_bcast(se_write_gll_corners, 1, mpi_logical, masterprocid, mpicom, ierr)
call MPI_bcast(se_horz_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
call MPI_bcast(se_vert_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
call MPI_bcast(se_tracer_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
call MPI_bcast(se_hypervis_on_plevs, 1, mpi_logical, masterprocid, mpicom, ierr)
call MPI_bcast(se_write_restart_unstruct, 1, mpi_logical, masterprocid, mpicom, ierr)
if (se_npes <= 0) then
call endrun('dyn_readnl: ERROR: se_npes must be > 0')
end if
! Initialize the SE structure that holds the MPI decomposition information
par = initmpi(se_npes)
call initomp()
!
! automatically set viscosity coefficients
!
!
! Use scaling from
!
! - Boville, B. A., 1991: Sensitivity of simulated climate to
! model resolution. J. Climate, 4, 469485.
!
! - TAKAHASHI ET AL., 2006: GLOBAL SIMULATION OF MESOSCALE SPECTRUM
!
uniform_res_hypervis_scaling = log(10.0_r8)/log(2.0_r8) !=3.3219
!
! compute factor so that at ne30 resolution nu=1E15
!
nu_fac = 1.0E15_r8/(110000.0_r8**uniform_res_hypervis_scaling)
if (se_nu_div < 0) then
if (se_ne <= 0) then
call endrun('dyn_readnl: ERROR must have se_ne > 0 for se_nu_div < 0')
end if
se_nu_div = nu_fac*((30.0_r8/se_ne)*110000.0_r8)**uniform_res_hypervis_scaling
end if
if (se_nu_p < 0) then
if (se_ne <= 0) then
call endrun('dyn_readnl: ERROR must have se_ne > 0 for se_nu_p < 0')
end if
se_nu_p = nu_fac*((30.0_r8/se_ne)*110000.0_r8)**uniform_res_hypervis_scaling
end if
if (se_nu < 0) then
if (se_ne <= 0) then
call endrun('dyn_readnl: ERROR must have se_ne > 0 for se_nu < 0')
end if
se_nu = 0.4_r8*nu_fac*((30.0_r8/se_ne)*110000.0_r8)**uniform_res_hypervis_scaling
end if
! Go ahead and enforce ne = 0 for refined mesh runs
if (se_refined_mesh) then
se_ne = 0
end if
! Set HOMME defaults
call homme_set_defaults()
! Set HOMME variables not in CAM's namelist but with different CAM defaults
partmethod = SFCURVE
npart = se_npes
! CAM requires forward-in-time, subcycled dynamics
! RK2 3 stage tracers, sign-preserving conservative
rk_stage_user = 3
topology = "cube"
! Finally, set the HOMME variables which have different names
qsize_condensate_loading = se_qsize_condensate_loading
lcp_moist = .false.!xxx
fine_ne = se_fine_ne
ftype = se_ftype
statediag_numtrac = MIN(se_statediag_numtrac,pcnst)
hypervis_power = se_hypervis_power
hypervis_scaling = se_hypervis_scaling
hypervis_subcycle = se_hypervis_subcycle
hypervis_subcycle_q = se_hypervis_subcycle_q
limiter_option = se_limiter_option
max_hypervis_courant = se_max_hypervis_courant
refined_mesh = se_refined_mesh
ne = se_ne
nsplit = se_nsplit
nu = se_nu
nu_div = se_nu_div
nu_p = se_nu_p
nu_q = se_nu_p !for tracer-wind consistency nu_q must me equal to nu_p
nu_top = se_nu_top
qsplit = se_qsplit
rsplit = se_rsplit
statefreq = se_statefreq
tstep_type = se_tstep_type
vert_remap_q_alg = se_vert_remap_q_alg
fv_nphys = se_fv_nphys
hypervis_on_plevs = se_hypervis_on_plevs
if (fv_nphys > 0) then
! Use finite volume physics grid and CSLAM for tracer advection
nphys_pts = fv_nphys*fv_nphys
tracer_transport_type = TRACERTRANSPORT_CONSISTENT_SE_FVM
qsize = qsize_condensate_loading ! number tracers advected by GLL
ntrac = pcnst ! number tracers advected by CSLAM
else
! Use GLL grid for physics and tracer advection
nphys_pts = npsq
tracer_transport_type = TRACERTRANSPORT_SE_GLL
qsize = pcnst
ntrac = 0
end if
if (rsplit < 1) then
call endrun('dyn_readnl: rsplit must be > 0')
end if
! if restart or branch run
if (.not. initial_run) then
runtype = 1
end if
! HOMME wants 'none' to indicate no mesh file
if (len_trim(se_mesh_file) == 0) then
se_mesh_file = 'none'
if (se_refined_mesh) then
call endrun('dyn_readnl ERROR: se_refined_mesh=.true. but no se_mesh_file')
end if
end if
call homme_postprocess_namelist(se_mesh_file, par)
! Set threading numbers to reasonable values
if ((se_horz_num_threads == 0) .and. (se_vert_num_threads == 0) .and. (se_tracer_num_threads == 0)) then
! The user has not set any threading values, choose defaults
se_horz_num_threads = max_num_threads
se_vert_num_threads = 1
se_tracer_num_threads = se_vert_num_threads
end if
if (se_horz_num_threads < 1) then
se_horz_num_threads = 1
end if
if (se_vert_num_threads < 1) then
se_vert_num_threads = 1
end if
if (se_tracer_num_threads < 1) then
se_tracer_num_threads = 1
end if
horz_num_threads = se_horz_num_threads
vert_num_threads = se_vert_num_threads
tracer_num_threads = se_tracer_num_threads
write_restart_unstruct = se_write_restart_unstruct
if (masterproc) then
write(iulog, '(a,i0)') 'dyn_readnl: se_ftype = ',ftype
write(iulog, '(a,i0)') 'dyn_readnl: se_statediag_numtrac = ',statediag_numtrac
write(iulog, '(a,i0)') 'dyn_readnl: se_hypervis_subcycle = ',se_hypervis_subcycle
write(iulog, '(a,i0)') 'dyn_readnl: se_hypervis_subcycle_q = ',se_hypervis_subcycle_q
write(iulog, '(a,i0)') 'dyn_readnl: se_limiter_option = ',se_limiter_option
if (.not. se_refined_mesh) then
write(iulog, '(a,i0)') 'dyn_readnl: se_ne = ',se_ne
end if
write(iulog, '(a,i0)') 'dyn_readnl: se_npes = ',se_npes
write(iulog, '(a,i0)') 'dyn_readnl: se_nsplit = ',se_nsplit
write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu = ',se_nu
write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_div = ',se_nu_div
write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_p = ',se_nu_p
write(iulog, '(a)') 'Note that nu_q=nu_p for mass / tracer inconsistency'
write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_top = ',se_nu_top
write(iulog, '(a,i0)') 'dyn_readnl: se_qsplit = ',se_qsplit
write(iulog, '(a,i0)') 'dyn_readnl: se_rsplit = ',se_rsplit
write(iulog, '(a,i0)') 'dyn_readnl: se_statefreq = ',se_statefreq
write(iulog, '(a,i0)') 'dyn_readnl: se_tstep_type = ',se_tstep_type
write(iulog, '(a,i0)') 'dyn_readnl: se_vert_remap_q_alg = ',se_vert_remap_q_alg
write(iulog, '(a,i0)') 'dyn_readnl: se_qsize_condensate_loading = ',se_qsize_condensate_loading
write(iulog, '(a,l4)') ' : lcp_moist = ',lcp_moist
write(iulog, '(a,l4)') ' dyn_readnl: hypervis_on_plevs = ',hypervis_on_plevs
if (hypervis_on_plevs.and.nu_p>0) then
write(iulog, *) 'FYI: nu_p>0 and hypervis_on_plevs=T => hypervis is applied to dp-dp_ref'
else if (hypervis_on_plevs.and.nu_p==0) then
write(iulog, *) 'FYI: hypervis_on_plevs=T and nu_p=0'
end if
if (se_refined_mesh) then
write(iulog, '(a)') 'dyn_readnl: Refined mesh simulation'
write(iulog, '(a)') 'dyn_readnl: se_mesh_file = ',trim(se_mesh_file)
if (abs(se_hypervis_power) < 1.0e-12_r8) then
write(iulog, '(a,e11.4)') 'dyn_readnl: se_hypervis_power = ',se_hypervis_power, ', (tensor hyperviscosity)'
write(iulog, '(a,e11.4)') 'dyn_readnl: se_hypervis_scaling = ',se_hypervis_scaling
else if (abs(se_hypervis_power - 3.322_r8) < 1.0e-12_r8) then
write(iulog, '(a,e11.4)') 'dyn_readnl: se_hypervis_power = ',se_hypervis_power, ', (scalar hyperviscosity)'
write(iulog, '(a,i0)') 'dyn_readnl: se_fine_ne = ',se_fine_ne
else
write(iulog, '(a,i0)') 'dyn_readnl: se_hypervis_power = ',se_hypervis_power
write(iulog, '(a,e11.4)') 'dyn_readnl: se_hypervis_scaling = ',se_hypervis_scaling
write(iulog, '(a,e11.4)') 'dyn_readnl: se_fine_ne = ',se_fine_ne
end if
write(iulog, '(a,e11.4)') 'dyn_readnl: se_max_hypervis_courant = ',se_max_hypervis_courant
end if
if ((se_met_nudge_u /= 0._r8) .or. (se_met_nudge_p /= 0._r8) .or. &
(se_met_nudge_t /= 0._r8) .or. (se_met_tevolve /= 0)) then
write(iulog, '(a)') 'dyn_readnl: Nudging:'
write(iulog,'(a,e14.6)') " : se_met_nudge_u = ", se_met_nudge_u
write(iulog,'(a,e14.6)') " : se_met_nudge_p = ", se_met_nudge_p
write(iulog,'(a,e14.6)') " : se_met_nudge_t = ", se_met_nudge_t
write(iulog,'(a,i0)') " : se_met_tevolve = ", se_met_tevolve
else
write(iulog, '(a)') 'dyn_readnl: Nudging off'
end if
if (fv_nphys > 0) then
write(iulog, '(a)') 'dyn_readnl: physics will run on FVM points; advection by CSLAM'
write(iulog,'(a,i0)') 'dyn_readnl: se_fv_nphys = ', fv_nphys
else
write(iulog, '(a)') 'dyn_readnl: physics will run on SE GLL points'
end if
write(iulog, '(a,i0)') 'dyn_readnl: se_horz_num_threads = ',horz_num_threads
write(iulog, '(a,i0)') 'dyn_readnl: se_vert_num_threads = ',vert_num_threads
write(iulog, '(a,i0)') 'dyn_readnl: se_tracer_num_threads = ',tracer_num_threads
if (trim(se_write_grid_file) == 'SCRIP') then
write(iulog,'(2a)') "dyn_readnl: write SCRIP grid file = ", trim(se_grid_filename)
else
write(iulog,'(a)') "dyn_readnl: do not write grid file"
end if
write(iulog,'(a,l1)') 'dyn_readnl: write gll corners to SEMapping.nc = ', &
se_write_gll_corners
write(iulog,'(a,l1)') 'dyn_readnl: write restart data on unstructured grid = ', &
se_write_restart_unstruct
end if
call native_mapping_readnl(NLFileName)
end subroutine dyn_readnl
!=========================================================================================
subroutine dyn_register()
use physics_buffer, only: pbuf_add_field, dtype_r8
use ppgrid, only: pcols, pver
! 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 dyn_grid, only: elem, fvm
use cam_pio_utils, only: clean_iodesc_list
use physconst, only: cpwv, cpliq, cpice
use cam_history, only: addfld, add_default, horiz_only, register_vector_field
use gravity_waves_sources, only: gws_init
use thread_mod, only: horz_num_threads
use hybrid_mod, only: get_loop_ranges, config_thread_region
use dimensions_mod, only: qsize_condensate_loading,qsize_condensate_loading_idx
use dimensions_mod, only: qsize_condensate_loading_idx_gll
use dimensions_mod, only: qsize_condensate_loading_cp
use dimensions_mod, only: cnst_name_gll, cnst_longname_gll
use prim_driver_mod, only: prim_init2
use time_mod, only: time_at
use control_mod, only: runtype
use test_fvm_mapping, only: test_mapping_addfld
! Dummy arguments:
type(dyn_import_t), intent(out) :: dyn_in
type(dyn_export_t), intent(out) :: dyn_out
! Local variables
integer :: ithr, nets, nete, ie, k
real(r8), parameter :: Tinit = 300.0_r8
type(hybrid_t) :: hybrid
integer :: ixcldice, ixcldliq, ixrain, ixsnow, ixgraupel
integer :: m_cnst, m
! variables for initializing energy and axial angular momentum diagnostics
character (len = 3), dimension(10) :: stage = (/"dED","dAF","dBD","dAD","dAR","dBF","dBH","dCH","dAH",'p2d'/)
character (len = 70),dimension(10) :: stage_txt = (/&
" end of previous dynamics ",& !dED
" from previous remapping or state passed to dynamics",& !dAF - state in beginning of nsplit loop
" state after applying CAM forcing ",& !dBD - state after applyCAMforcing
" before vertical remapping ",& !dAD - state before vertical remapping
" after vertical remapping ",& !dAR - state at end of nsplit loop
" state passed to parameterizations ",& !dBF
" state before hypervis ",& !dBH
" state after hypervis but before adding heating term",& !dCH
" state after hypervis ",& !dAH
" phys2dyn mapping errors (requires ftype-1) " & !p2d - for assessing phys2dyn mapping errors
/)
character (len = 2) , dimension(8) :: vars = (/"WV","WL","WI","SE","KE","MR","MO","TT"/)
character (len = 70), dimension(8) :: vars_descriptor = (/&
"Total column water vapor ",&
"Total column cloud water ",&
"Total column cloud ice ",&
"Total column dry static energy ",&
"Total column kinetic energy ",&
"Total column wind axial angular momentum",&
"Total column mass axial angular momentum",&
"Total column test tracer "/)
character (len = 14), dimension(8) :: &
vars_unit = (/&
"kg/m2 ","kg/m2 ","kg/m2 ","J/m2 ",&
"J/m2 ","kg*m2/s*rad2 ","kg*m2/s*rad2 ","kg/m2 "/)
integer :: istage, ivars
character (len=108) :: str1, str2, str3
integer, parameter :: qcondensate_max = 6
character(len=*), parameter :: subname = 'dyn_init'
!----------------------------------------------------------------------------
if (qsize_condensate_loading > qcondensate_max) then
call endrun(subname//': se_qsize_condensate_loading not setup for more than 6 forms of water')
end if
! Now allocate and set condenstate vars
allocate(qsize_condensate_loading_idx(qcondensate_max))
allocate(qsize_condensate_loading_idx_gll(qcondensate_max))
allocate(qsize_condensate_loading_cp(qcondensate_max))
allocate(cnst_name_gll(qsize)) ! constituent names for gll tracers
allocate(cnst_longname_gll(qsize)) ! long name of constituents for gll tracers
! water vapor is always tracer 1
qsize_condensate_loading_idx(1) = 1
qsize_condensate_loading_cp(1) = cpwv
call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
if (ixcldliq < 1.and.qsize_condensate_loading > 1) &
call endrun(subname//': ERROR: qsize_condensate_loading >1 but CLDLIQ not available')
qsize_condensate_loading_idx(2) = ixcldliq
qsize_condensate_loading_cp(2) = cpliq
call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
if (ixcldice < 1.and.qsize_condensate_loading > 2) &
call endrun(subname//': ERROR: qsize_condensate_loading >2 but CLDICE not available')
qsize_condensate_loading_idx(3) = ixcldice
qsize_condensate_loading_cp(3) = cpice
call cnst_get_ind('RAINQM', ixrain, abort=.false.)
if (ixrain < 1.and.qsize_condensate_loading > 3) &
call endrun(subname//': ERROR: qsize_condensate_loading >3 but RAINQM not available')
qsize_condensate_loading_idx(4) = ixrain
qsize_condensate_loading_cp(4) = cpliq
call cnst_get_ind('SNOWQM', ixsnow, abort=.false.)
if (ixsnow < 1.and.qsize_condensate_loading > 4) &
call endrun(subname//': ERROR: qsize_condensate_loading >4 but SNOWQM not available')
qsize_condensate_loading_idx(5) = ixsnow
qsize_condensate_loading_cp(5) = cpice
call cnst_get_ind('GRAUQM', ixgraupel, abort=.false.)
if (ixgraupel < 1.and.qsize_condensate_loading > 5) &
call endrun(subname//': ERROR: qsize_condensate_loading >5 but GRAUQM not available')
qsize_condensate_loading_idx(6) = ixgraupel
qsize_condensate_loading_cp(6) = cpice
!
! if adding more condensate loading tracers remember to increase qsize_d in dimensions_mod
!
qsize_condensate_loading_idx_gll(:) = -1
do m=1,qsize
!
! The "_gll" index variables below are used to keep track of condensate-loading tracers
! since they are not necessarily indexed contiguously and not necessarily in the same
! order (physics is in charge of the order)
!
! if running with CSLAM then the SE (gll) condensate-loading water tracers are always
! indexed contiguously (q,cldliq,cldice,rain,snow,graupel) - see above
!
! CSLAM tracers are always indexed as in physics
! of no CSLAM then SE tracers are always indexed as in physics
!
if (ntrac>0) then
!
! note that in this case qsize = qsize_condensate_loading
!
qsize_condensate_loading_idx_gll(m) = m
cnst_name_gll (m) = cnst_name (qsize_condensate_loading_idx(m))
cnst_longname_gll(m) = cnst_longname(qsize_condensate_loading_idx(m))
else
!
! if not running with CSLAM then the condensate-loading water tracers are not necessarily
! indexed contiguously (are indexed as in physics)
!
if (m.le.qcondensate_max) qsize_condensate_loading_idx_gll(m) = qsize_condensate_loading_idx(m)
cnst_name_gll (m) = cnst_name (m)
cnst_longname_gll(m) = cnst_longname(m)
end if
end do
! if user wants to add more condensate loading tracers add them here ....
!
! Initialize the import/export objects
!
if(iam < par%nprocs) then
dyn_in%elem => elem
dyn_in%fvm => fvm
dyn_out%elem => elem
dyn_out%fvm => fvm
else
nullify(dyn_in%elem)
nullify(dyn_in%fvm)
nullify(dyn_out%elem)
nullify(dyn_out%fvm)
end if
call read_phis(dyn_in)
if (initial_run) then
call read_inidat(dyn_in)
call clean_iodesc_list()
end if
if (iam < par%nprocs) then
!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,ie)
hybrid = config_thread_region(par,'horizontal')
call get_loop_ranges(hybrid, ibeg=nets, iend=nete)
call prim_init2(elem, fvm, hybrid, nets, nete, TimeLevel, hvcoord)
!$OMP END PARALLEL
if (use_gw_front .or. use_gw_front_igw) call gws_init(elem)
end if ! iam < par%nprocs
! Forcing from physics on the GLL grid
call addfld ('FU', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind forcing term on GLL grid', gridname='GLL')
call addfld ('FV', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind forcing term on GLL grid',gridname='GLL')
call register_vector_field('FU', 'FV')
call addfld ('FT', (/ 'lev' /), 'A', 'K/s', 'Temperature forcing term on GLL grid',gridname='GLL')
! Tracer forcing on fvm (CSLAM) grid and internal CSLAM pressure fields
if (ntrac>0) then
do m = 1, ntrac
call addfld (trim(cnst_name(m))//'_fvm', (/ 'lev' /), 'I', 'kg/kg', &
trim(cnst_longname(m)), gridname='FVM')
call addfld ('F'//trim(cnst_name(m))//'_fvm', (/ 'lev' /), 'I', 'kg/kg/s', &
trim(cnst_longname(m))//' mixing ratio forcing term (q_new-q_old) on fvm grid', &
gridname='FVM')
end do
call addfld ('dp_fvm' ,(/ 'lev' /), 'I', 'Pa','CSLAM Pressure level thickness', gridname='FVM')
call addfld ('PSDRY_fvm',horiz_only, 'I','Pa','CSLAM dry surface pressure' , gridname='FVM')
end if
do m_cnst = 1, qsize
call addfld ('F'//trim(cnst_name_gll(m_cnst))//'_gll', (/ 'lev' /), 'I', 'kg/kg/s', &
trim(cnst_longname(m_cnst))//' mixing ratio forcing term (q_new-q_old) on GLL grid', gridname='GLL')
end do
! Energy diagnostics and axial angular momentum diagnostics
call addfld ('ABS_dPSdt', horiz_only, 'A', 'Pa/s', 'Absolute surface pressure tendency',gridname='GLL')
if (ntrac>0) then
call addfld ('WV_PDC', horiz_only, 'A', 'kg/m2','Total column water vapor lost in physics-dynamics coupling',gridname='FVM')
call addfld ('WL_PDC', horiz_only, 'A', 'kg/m2','Total column cloud water lost in physics-dynamics coupling',gridname='FVM')
call addfld ('WI_PDC', horiz_only, 'A', 'kg/m2','Total column cloud ice lost in physics-dynamics coupling' ,gridname='FVM')
call addfld ('TT_PDC', horiz_only, 'A', 'kg/m2','Total column test tracer lost in physics-dynamics coupling',gridname='FVM')
else
call addfld ('WV_PDC', horiz_only, 'A', 'kg/m2','Total column water vapor lost in physics-dynamics coupling',gridname='GLL')
call addfld ('WL_PDC', horiz_only, 'A', 'kg/m2','Total column cloud water lost in physics-dynamics coupling',gridname='GLL')
call addfld ('WI_PDC', horiz_only, 'A', 'kg/m2','Total column cloud ice lost in physics-dynamics coupling' ,gridname='GLL')
call addfld ('TT_PDC', horiz_only, 'A', 'kg/m2','Total column test tracer lost in physics-dynamics coupling',gridname='GLL')
end if
do istage = 1,SIZE(stage)
do ivars=1,SIZE(vars)
write(str1,*) TRIM(ADJUSTL(vars(ivars))),TRIM(ADJUSTL("_")),TRIM(ADJUSTL(stage(istage)))
write(str2,*) TRIM(ADJUSTL(vars_descriptor(ivars))),&
TRIM(ADJUSTL(" ")),TRIM(ADJUSTL(stage_txt(istage)))
write(str3,*) TRIM(ADJUSTL(vars_unit(ivars)))
call addfld (TRIM(ADJUSTL(str1)), horiz_only, 'A', TRIM(ADJUSTL(str3)),TRIM(ADJUSTL(str2)),gridname='GLL')
end do
end do
call test_mapping_addfld
#ifdef extra_pdc_diags
!xxx
!xxx this is just for debugging - remove for trunk
!xxx
if (ntrac>0) then
call addfld ('WV_physgrid', (/ 'lev' /), 'A', 'kg/m2','Total column water vapor forcing physgrid')
call addfld ('WV_fvm', (/ 'lev' /), 'A', 'kg/m2','Total column water vapor forcing fvm',gridname='FVM')
call addfld ('WL_physgrid', (/ 'lev' /), 'A', 'kg/m2','Total column cldliq forcing physgrid')
call addfld ('WL_fvm', (/ 'lev' /), 'A', 'kg/m2','Total column cldliq forcing fvm',gridname='FVM')
call addfld ('WI_physgrid', (/ 'lev' /), 'A', 'kg/m2','Total column cldice forcing physgrid')
call addfld ('WI_fvm', (/ 'lev' /), 'A', 'kg/m2','Total column cldice forcing fvm',gridname='FVM')
call addfld ('TT_physgrid', (/ 'lev' /), 'A', 'kg/m2','Total column TT forcing physgrid')
call addfld ('TT_fvm', (/ 'lev' /), 'A', 'kg/m2','Total column TT forcing fvm',gridname='FVM')
end if
#endif
end subroutine dyn_init
!=========================================================================================
subroutine dyn_run(dyn_state)
use prim_advance_mod, only: calc_tot_energy_dynamics
use prim_driver_mod, only: prim_run_subcycle
use dimensions_mod, only: qsize_condensate_loading, qsize_condensate_loading_idx_gll
use dimensions_mod, only: cnst_name_gll
use time_mod, only: tstep, nsplit, timelevel_qdp
use hybrid_mod, only: config_thread_region, get_loop_ranges
use control_mod, only: qsplit ,rsplit
use thread_mod, only: horz_num_threads
use time_mod, only: tevolve
type(dyn_export_t), intent(inout) :: dyn_state
type(hybrid_t) :: hybrid
integer :: tl_f
integer :: n
integer :: nets, nete, ithr
integer :: i, ie, j, k, m
integer :: n0_qdp
logical :: ldiag
real(r8) :: ftmp(npsq,nlev,3)
real(r8) :: dtime
real(r8) :: rec2dt
real(r8), allocatable, dimension(:,:,:) :: ps_before
real(r8), allocatable, dimension(:,:,:) :: abs_ps_tend
!----------------------------------------------------------------------------
#ifdef debug_coupling
return
#endif
tevolve = 0._r8
if (iam >= par%nprocs) return
!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n,ie,m,i,j,k)
hybrid = config_thread_region(par,'horizontal')
call get_loop_ranges(hybrid, ibeg=nets, iend=nete)
dtime = get_step_size()
rec2dt = 1._r8/dtime
tl_f = TimeLevel%n0 ! timelevel which was adjusted by physics
! output physics forcing
if (hist_fld_active('FU') .or. hist_fld_active('FV') .or.hist_fld_active('FT')) then
do ie = nets, nete
do k = 1, nlev
do j = 1, np
do i = 1, np
ftmp(i+(j-1)*np,k,1) = dyn_state%elem(ie)%derived%FM(i,j,1,k)
ftmp(i+(j-1)*np,k,2) = dyn_state%elem(ie)%derived%FM(i,j,2,k)
ftmp(i+(j-1)*np,k,3) = dyn_state%elem(ie)%derived%FT(i,j,k)
end do
end do
end do
call outfld('FU', ftmp(:,:,1), npsq, ie)
call outfld('FV', ftmp(:,:,2), npsq, ie)
call outfld('FT', ftmp(:,:,3), npsq, ie)
end do
end if
do m = 1, qsize
if (hist_fld_active('F'//trim(cnst_name_gll(m))//'_gll')) then
do ie = nets, nete
call outfld('F'//trim(cnst_name_gll(m))//'_gll',&
RESHAPE(dyn_state%elem(ie)%derived%FQ(:,:,:,m), (/np*np,nlev/)), npsq, ie)
end do
end if
end do
! convert elem(ie)%derived%fq to tendency
do ie = nets, nete
do m = 1, qsize
do k = 1, nlev
do j = 1, np
do i = 1, np
dyn_state%elem(ie)%derived%FQ(i,j,k,m) = dyn_state%elem(ie)%derived%FQ(i,j,k,m)* &
rec2dt*dyn_state%elem(ie)%state%dp3d(i,j,k,tl_f)
end do
end do
end do
end do
end do
if (ntrac > 0) then
do ie = nets, nete
do m = 1, ntrac
do k = 1, nlev
do j = 1, nc
do i = 1, nc
dyn_state%fvm(ie)%fc(i,j,k,m) = dyn_state%fvm(ie)%fc(i,j,k,m)* &
rec2dt!*dyn_state%fvm(ie)%dp_fvm(i,j,k,n0_fvm)
end do
end do
end do
end do
end do
end if
ldiag = hist_fld_active('ABS_dPSdt')
if (ldiag) then
allocate(ps_before(np,np,nets:nete))
allocate(abs_ps_tend(np,np,nets:nete))
abs_ps_tend(:,:,nets:nete) = 0.0_r8
end if
do n = 1, nsplit
if (ldiag) then
do ie = nets, nete
ps_before(:,:,ie) = dyn_state%elem(ie)%state%psdry(:,:,tl_f)
end do
end if
! forward-in-time RK, with subcycling
call prim_run_subcycle(dyn_state%elem, dyn_state%fvm, hybrid, nets, nete, &
tstep, TimeLevel, hvcoord, n)
if (ldiag) then
do ie = nets, nete
abs_ps_tend(:,:,ie) = abs_ps_tend(:,:,ie) + &
ABS(ps_before(:,:,ie)-dyn_state%elem(ie)%state%psdry(:,:,tl_f)) &
/(tstep*qsplit*rsplit)
end do
end if
end do
if (ldiag) then
do ie=nets,nete
abs_ps_tend(:,:,ie)=abs_ps_tend(:,:,ie)/DBLE(nsplit)
call outfld('ABS_dPSdt',RESHAPE(abs_ps_tend(:,:,ie),(/npsq/)),npsq,ie)
end do
deallocate(ps_before,abs_ps_tend)
end if
call TimeLevel_Qdp(TimeLevel, qsplit, n0_qdp)!get n0_qdp for diagnostics call
call calc_tot_energy_dynamics(dyn_state%elem, nets, nete, tl_f, n0_qdp, 'dBF')
!$OMP END PARALLEL
! output vars on CSLAM fvm grid
call write_dyn_vars(dyn_state)
end subroutine dyn_run
!===============================================================================
subroutine dyn_final(DYN_STATE, RESTART_FILE)
type (elem_state_t), target :: DYN_STATE
character(LEN=*) , intent(IN) :: RESTART_FILE
end subroutine dyn_final
!===============================================================================
subroutine read_inidat(dyn_in)
use shr_sys_mod, only: shr_sys_flush
use hycoef, only: hyai, hybi, ps0
use const_init, only: cnst_init_default
use cam_control_mod, only: ideal_phys
use element_mod, only: timelevels
use dimensions_mod, only: qsize_d, qsize_condensate_loading
use dimensions_mod, only: qsize_condensate_loading_idx
use fvm_mapping, only: dyn2fvm_mass_vars
use control_mod, only: runtype,initial_global_ave_dry_ps
use prim_driver_mod, only: prim_set_dry_mass
! Arguments
type (dyn_import_t), target, intent(inout) :: dyn_in ! dynamics import
! Local variables
integer(iMap), pointer :: ldof(:) ! Basic (2D) grid dof
type(file_desc_t), pointer :: fh_ini, fh_topo
type(element_t), pointer :: elem(:)
real(r8), allocatable :: qtmp(:,:,:,:,:) ! (np,np,nlev,nelemd,n)
real(r8), allocatable :: dbuf2(:,:) ! (npsq,nelemd)
real(r8), allocatable :: dbuf3(:,:,:) ! (npsq,nlev,nelemd)
real(r8), allocatable :: factor_array(:,:,:,:) ! (np,np,nlev,nelemd)
logical, allocatable :: pmask(:) ! (npsq*nelemd) unique grid vals
logical, allocatable :: pmask_phys(:)
character(len=max_hcoordname_len):: grid_name
real(r8), allocatable :: latvals(:),latvals_phys(:)
real(r8), allocatable :: lonvals(:),lonvals_phys(:)
real(r8), pointer :: latvals_deg(:)
real(r8), pointer :: lonvals_deg(:)
integer :: ie, k, t
character(len=max_fieldname_len) :: fieldname, fieldname2
logical :: found
logical :: inic_wet ! true if initial condition is based on
! wet pressure and water species
integer :: kptr, m_cnst
type(EdgeBuffer_t) :: edge
integer :: lsize
character(len=max_fieldname_len) :: dimname, varname
integer :: ierr
integer :: ncol_did
integer :: ncol_size
integer :: rndm_seed_sz
integer, allocatable :: rndm_seed(:)
integer :: dims(2)
integer :: pio_errtype
real(r8) :: pertval
integer :: i, j, indx, nq
integer :: dyn_cols
character(len=128) :: errmsg
character(len=*), parameter :: subname='READ_INIDAT'
integer :: ioff
! fvm vars
real(r8), allocatable :: inv_dp_darea_fvm(:,:,:)
real(r8) :: min_val, max_val
real(r8) :: dp_tmp, pstmp(np,np)
! Variables for analytic initial conditions
integer, allocatable :: glob_ind(:)
integer, allocatable :: m_ind(:)
real(r8), allocatable :: dbuf4(:,:,:,:)
!----------------------------------------------------------------------------
fh_ini => initial_file_get_id()
fh_topo => topo_file_get_id()