-
Notifications
You must be signed in to change notification settings - Fork 12
/
grid_class.F90
2339 lines (1873 loc) · 81.9 KB
/
grid_class.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
!------------------------------------------------------------------------
! $Id$
!===============================================================================
module grid_class
! Description:
!
! Definition of a grid class and associated functions
!
! The grid specification is as follows:
!
! + ================== zm(nz) =========GP=========
! |
! |
! 1/dzt(nz) + ------------------ zt(nz) ---------GP---------
! | |
! | |
! + 1/dzm(nz-1) ================== zm(nz-1) ==================
! |
! |
! + ------------------ zt(nz-1) ------------------
!
! .
! .
! .
! .
!
! ================== zm(k+1) ===================
!
!
! + ------------------ zt(k+1) -------------------
! |
! |
! + 1/dzm(k) ================== zm(k) =====================
! | |
! | |
! 1/dzt(k) + ------------------ zt(k) ---------------------
! |
! |
! + ================== zm(k-1) ===================
!
!
! ------------------ zt(k-1) -------------------
!
! .
! .
! .
! .
!
! + ================== zm(2) =====================
! |
! |
! 1/dzt(2) + ------------------ zt(2) ---------------------
! | |
! | |
! + 1/dzm(1) ================== zm(1) ============GP======= zm_init
! | ////////////////////////////////////////////// surface
! |
! + ------------------ zt(1) ------------GP-------
!
!
! The variable zm(k) stands for the momentum level altitude at momentum
! level k; the variable zt(k) stands for the thermodynamic level altitude at
! thermodynamic level k; the variable invrs_dzt(k) is the inverse distance
! between momentum levels (over a central thermodynamic level k); and the
! variable invrs_dzm(k) is the inverse distance between thermodynamic levels
! (over a central momentum level k). Please note that in the above diagram,
! "invrs_dzt" is denoted "dzt", and "invrs_dzm" is denoted "dzm", such that
! 1/dzt is the distance between successive momentum levels k-1 and k (over a
! central thermodynamic level k), and 1/dzm is the distance between successive
! thermodynamic levels k and k+1 (over a central momentum level k).
!
! The grid setup is compatible with a stretched (unevely-spaced) grid. Thus,
! the distance between successive grid levels may not always be constant.
!
! The following diagram is an example of a stretched grid that is defined on
! momentum levels. The thermodynamic levels are placed exactly halfway
! between the momentum levels. However, the momentum levels do not fall
! halfway between the thermodynamic levels.
!
! =============== zm(k+1) ===============
!
!
!
! --------------- zt(k+1) ---------------
!
!
!
! =============== zm(k) ===============
!
! --------------- zt(k) ---------------
!
! =============== zm(k-1) ===============
!
! The following diagram is an example of a stretched grid that is defined on
! thermodynamic levels. The momentum levels are placed exactly halfway
! between the thermodynamic levels. However, the thermodynamic levels do not
! fall halfway between the momentum levels.
!
! --------------- zt(k+1) ---------------
!
!
!
! =============== zm(k) ===============
!
!
!
! --------------- zt(k) ---------------
!
! =============== zm(k-1) ===============
!
! --------------- zt(k-1) ---------------
!
! NOTE: Any future code written for use in the CLUBB parameterization should
! use interpolation formulas consistent with a stretched grid. The
! simplest way to do so is to call the appropriate interpolation
! function from this module. Interpolations should *not* be handled in
! the form of: ( var_zm(k) + var_zm(k-1) ) / 2; *nor* in the form of:
! 0.5_core_rknd*( var_zt(k+1) + var_zt(k) ). Rather, all explicit interpolations
! should call zt2zm or zm2zt; while interpolations for a variable being
! solved for implicitly in the code should use gr%weights_zt2zm (which
! refers to interp_weights_zt2zm_imp), or gr%weights_zm2zt (which
! refers to interp_weights_zm2zt_imp).
!
! Momentum level 1 is placed at altitude zm_init, which is usually at the
! surface. However, in general, zm_init can be at any altitude defined by the
! user.
!
! GP indicates ghost points. Variables located at those levels are not
! prognosed, but only used for boundary conditions.
!
! Chris Golaz, 7/17/99
! modified 9/10/99
! schemena, modified 6/11/2014 - Restructered code to add cubic/linear flag
! References:
!
! https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:clubb_grid
!
! Section 3c, p. 3548 /Numerical discretization/ of:
! ``A PDF-Based Model for Boundary Layer Clouds. Part I:
! Method and Model Description'' Golaz, et al. (2002)
! JAS, Vol. 59, pp. 3540--3551.
!-----------------------------------------------------------------------
use clubb_precision, only: &
core_rknd ! Variable(s)
implicit none
public :: gr, grid, zt2zm, interp_weights_zt2zm_imp, zm2zt, &
interp_weights_zm2zt_imp, ddzm, ddzt, &
setup_grid, cleanup_grid, setup_grid_heights, &
read_grid_heights, flip
private :: linear_interpolated_azm, linear_interpolated_azmk, &
interpolated_azmk_imp, linear_interpolated_azt, &
linear_interpolated_aztk, interpolated_aztk_imp, &
gradzm, gradzt, t_above, t_below, m_above, m_below, &
cubic_interpolated_azmk, cubic_interpolated_aztk, &
cubic_interpolated_azm, cubic_interpolated_azt
private ! Default Scoping
! Constant parameters
integer, parameter :: &
t_above = 1, & ! Upper thermodynamic level index (gr%weights_zt2zm).
t_below = 2, & ! Lower thermodynamic level index (gr%weights_zt2zm).
m_above = 1, & ! Upper momentum level index (gr%weights_zm2zt).
m_below = 2 ! Lower momentum level index (gr%weights_zm2zt).
type grid
integer :: nz ! Number of points in the grid
! Note: Fortran 90/95 prevents an allocatable array from appearing
! within a derived type. However, Fortran 2003 does not!!!!!!!!
real( kind = core_rknd ), allocatable, dimension(:) :: &
zm, & ! Momentum grid
zt ! Thermo grid
real( kind = core_rknd ), allocatable, dimension(:) :: &
invrs_dzm, & ! The inverse spacing between thermodynamic grid
! levels; centered over momentum grid levels.
invrs_dzt ! The inverse spacing between momentum grid levels;
! centered over thermodynamic grid levels.
real( kind = core_rknd ), allocatable, dimension(:) :: &
dzm, & ! Spacing between thermodynamic grid levels; centered over
! momentum grid levels
dzt ! Spcaing between momentum grid levels; centered over
! thermodynamic grid levels
! These weights are normally used in situations
! where a momentum level variable is being
! solved for implicitly in an equation and
! needs to be interpolated to the thermodynamic grid levels.
real( kind = core_rknd ), allocatable, dimension(:,:) :: weights_zm2zt, &
! These weights are normally used in situations where a
! thermodynamic level variable is being solved for implicitly in an equation
! and needs to be interpolated to the momentum grid levels.
weights_zt2zm
end type grid
! The grid is defined here so that it is common throughout the module.
! The implication is that only one grid can be defined !
type (grid), target :: gr
! Modification for using CLUBB in a host model (i.e. one grid per column)
!$omp threadprivate(gr)
! Interfaces provided for function overloading
interface zt2zm
! For l_cubic_interp = .true.
! This version uses cublic spline interpolation of Stefen (1990).
!
! For l_cubic_interp = .false.
! This performs a linear extension at the highest grid level and therefore
! does not guarantee, for positive definite quantities (e.g. wp2), that the
! extended point is indeed positive definite. Positive definiteness can be
! ensured with a max statement.
! In the future, we could add a flag (lposdef) and, when needed, apply the
! max statement directly within interpolated_azm and interpolated_azmk.
module procedure redirect_interpolated_azmk, redirect_interpolated_azm
end interface
interface zm2zt
! For l_cubic_interp = .true.
! This version uses cublic spline interpolation of Stefen (1990).
!
! For l_cubic_interp = .false.
! This performs a linear extension at the lowest grid level and therefore
! does not guarantee, for positive definite quantities (e.g. wp2), that the
! extended point is indeed positive definite. Positive definiteness can be
! ensured with a max statement.
! In the future, we could add a flag (lposdef) and, when needed, apply the
! max statement directly within interpolated_azt and interpolated_aztk.
module procedure redirect_interpolated_aztk, redirect_interpolated_azt
end interface
interface interp_weights_zt2zm_imp
module procedure interpolated_azmk_imp
end interface
interface interp_weights_zm2zt_imp
module procedure interpolated_aztk_imp
end interface
! Vertical derivative functions
interface ddzm
module procedure gradzm
end interface
interface ddzt
module procedure gradzt
end interface
contains
!=============================================================================
subroutine setup_grid( nzmax, sfc_elevation, l_implemented, &
grid_type, deltaz, zm_init, zm_top, &
momentum_heights, thermodynamic_heights, &
begin_height, end_height )
! Description:
! Grid Constructor
!
! This subroutine sets up the CLUBB vertical grid.
!
! References:
! ``Equations for CLUBB'', Sec. 8, Grid Configuration.
!-----------------------------------------------------------------------
use constants_clubb, only: &
fstderr ! Variable(s)
use error_code, only: &
clubb_at_least_debug_level, & ! Procedure
err_code, & ! Error indicator
clubb_fatal_error, & ! Constant
err_header ! String
use clubb_precision, only: &
core_rknd ! Variable(s)
implicit none
! Constant parameters
integer, parameter :: &
NWARNING = 250 ! Issue a warning if nzmax exceeds this number.
! Input Variables
integer, intent(in) :: &
nzmax ! Number of vertical levels in grid [#]
real( kind = core_rknd ), intent(in) :: &
sfc_elevation ! Elevation of ground level [m AMSL]
! Flag to see if CLUBB is running on it's own,
! or if it's implemented as part of a host model.
logical, intent(in) :: l_implemented
! If CLUBB is running on it's own, this option determines if it is using:
! 1) an evenly-spaced grid;
! 2) a stretched (unevenly-spaced) grid entered on the thermodynamic grid
! levels (with momentum levels set halfway between thermodynamic levels);
! or
! 3) a stretched (unevenly-spaced) grid entered on the momentum grid levels
! (with thermodynamic levels set halfway between momentum levels).
integer, intent(in) :: grid_type
! If the CLUBB model is running by itself, and is using an evenly-spaced
! grid (grid_type = 1), it needs the vertical grid spacing and
! momentum-level starting altitude as input.
real( kind = core_rknd ), intent(in) :: &
deltaz, & ! Vertical grid spacing [m]
zm_init, & ! Initial grid altitude (momentum level) [m]
zm_top ! Maximum grid altitude (momentum level) [m]
! If the CLUBB parameterization is implemented in a host model, it needs to
! use the host model's momentum level altitudes and thermodynamic level
! altitudes.
! If the CLUBB model is running by itself, but is using a stretched grid
! entered on thermodynamic levels (grid_type = 2), it needs to use the
! thermodynamic level altitudes as input.
! If the CLUBB model is running by itself, but is using a stretched grid
! entered on momentum levels (grid_type = 3), it needs to use the momentum
! level altitudes as input.
real( kind = core_rknd ), intent(in), dimension(nzmax) :: &
momentum_heights, & ! Momentum level altitudes (input) [m]
thermodynamic_heights ! Thermodynamic level altitudes (input) [m]
integer, intent(out) :: &
begin_height, & ! Lower bound for *_heights arrays [-]
end_height ! Upper bound for *_heights arrays [-]
! Local Variables
integer :: ierr, & ! Allocation stat
i ! Loop index
! ---- Begin Code ----
! Define the grid size
if ( nzmax > NWARNING .and. clubb_at_least_debug_level( 1 ) ) then
write(fstderr,*) "Warning: running with vertical grid "// &
"which is larger than", NWARNING, "levels."
write(fstderr,*) "This may take a lot of CPU time and memory."
end if
gr%nz = nzmax
! Default bounds
begin_height = 1
end_height = gr%nz
!---------------------------------------------------
if ( .not. l_implemented ) then
if ( grid_type == 1 ) then
! Determine the number of grid points given the spacing
! to fit within the bounds without going over.
gr%nz = floor( ( zm_top - zm_init + deltaz ) / deltaz )
else if( grid_type == 2 ) then! Thermo
! Find begin_height (lower bound)
i = gr%nz
do while( thermodynamic_heights(i) >= zm_init .and. i > 1 )
i = i - 1
end do
if( thermodynamic_heights(i) >= zm_init ) then
write(fstderr,*) err_header, "Stretched zt grid cannot fulfill zm_init requirement"
err_code = clubb_fatal_error
return
else
begin_height = i
end if
! Find end_height (upper bound)
i = gr%nz
do while( thermodynamic_heights(i) > zm_top .and. i > 1 )
i = i - 1
end do
if( zm_top < thermodynamic_heights(i) ) then
write(fstderr,*) err_header, "Stretched zt grid cannot fulfill zm_top requirement"
err_code = clubb_fatal_error
return
else
end_height = i
gr%nz = size( thermodynamic_heights(begin_height:end_height) )
end if
else if( grid_type == 3 ) then ! Momentum
! Find begin_height (lower bound)
i = 1
do while( momentum_heights(i) < zm_init .and. i < gr%nz )
i = i + 1
end do
if( momentum_heights(i) < zm_init ) then
write(fstderr,*) err_header, "Stretched zm grid cannot fulfill zm_init requirement"
err_code = clubb_fatal_error
return
else
begin_height = i
end if
! Find end_height (lower bound)
i = gr%nz
do while( momentum_heights(i) > zm_top .and. i > 1 )
i = i - 1
end do
if( momentum_heights(i) > zm_top ) then
write(fstderr,*) err_header, "Stretched zm grid cannot fulfill zm_top requirement"
err_code = clubb_fatal_error
return
else
end_height = i
gr%nz = size( momentum_heights(begin_height:end_height) )
end if
endif ! grid_type
endif ! .not. l_implemented
!---------------------------------------------------
! Allocate memory for the grid levels
allocate( gr%zm(gr%nz), gr%zt(gr%nz), &
gr%dzm(gr%nz), gr%dzt(gr%nz), &
gr%invrs_dzm(gr%nz), gr%invrs_dzt(gr%nz), &
gr%weights_zm2zt(m_above:m_below,gr%nz), &
gr%weights_zt2zm(t_above:t_below,gr%nz), &
stat=ierr )
if ( ierr /= 0 ) then
write(fstderr,*) err_header, "In setup_grid: allocation of grid variables failed."
err_code = clubb_fatal_error
return
end if
! Set the values for the derived types used for heights, derivatives, and
! interpolation from the momentum/thermodynamic grid
call setup_grid_heights &
( l_implemented, grid_type, &
deltaz, zm_init, &
momentum_heights(begin_height:end_height), &
thermodynamic_heights(begin_height:end_height) )
if ( sfc_elevation > gr%zm(1) ) then
write(fstderr,*) "The altitude of the lowest momentum level, " &
// "gr%zm(1), must be at or above the altitude of " &
// "the surface, sfc_elevation. The lowest model " &
// "momentum level cannot be below the surface."
write(fstderr,*) "Altitude of lowest momentum level =", gr%zm(1)
write(fstderr,*) "Altitude of the surface =", sfc_elevation
err_code = clubb_fatal_error
return
endif
return
end subroutine setup_grid
!=============================================================================
subroutine cleanup_grid
! Description:
! De-allocates the memory for the grid
!
! References:
! None
!------------------------------------------------------------------------------
use constants_clubb, only: &
fstderr ! Constant(s)
implicit none
! Local Variable(s)
integer :: ierr
! ----- Begin Code -----
! Allocate memory for grid levels
deallocate( gr%zm, gr%zt, &
gr%dzm, gr%dzt, &
gr%invrs_dzm, gr%invrs_dzt, &
gr%weights_zm2zt, gr%weights_zt2zm, &
stat=ierr )
if ( ierr /= 0 ) then
write(fstderr,*) "Grid deallocation failed."
end if
return
end subroutine cleanup_grid
!=============================================================================
subroutine setup_grid_heights &
( l_implemented, grid_type, &
deltaz, zm_init, momentum_heights, &
thermodynamic_heights )
! Description:
! Sets the heights and interpolation weights of the column.
! This is seperated from setup_grid for those host models that have heights
! that vary with time.
! References:
! None
!------------------------------------------------------------------------------
use constants_clubb, only: &
fstderr ! Constant(s)
use clubb_precision, only: &
core_rknd ! Variable(s)
use error_code, only: &
err_code, & ! Error indicator
clubb_fatal_error, & ! Constant
err_header ! String
implicit none
! Input Variables
! Flag to see if CLUBB is running on it's own,
! or if it's implemented as part of a host model.
logical, intent(in) :: l_implemented
! If CLUBB is running on it's own, this option determines if it is using:
! 1) an evenly-spaced grid;
! 2) a stretched (unevenly-spaced) grid entered on the thermodynamic grid
! levels (with momentum levels set halfway between thermodynamic levels);
! or
! 3) a stretched (unevenly-spaced) grid entered on the momentum grid levels
! (with thermodynamic levels set halfway between momentum levels).
integer, intent(in) :: grid_type
! If the CLUBB model is running by itself, and is using an evenly-spaced
! grid (grid_type = 1), it needs the vertical grid spacing and
! momentum-level starting altitude as input.
real( kind = core_rknd ), intent(in) :: &
deltaz, & ! Vertical grid spacing [m]
zm_init ! Initial grid altitude (momentum level) [m]
! If the CLUBB parameterization is implemented in a host model, it needs to
! use the host model's momentum level altitudes and thermodynamic level
! altitudes.
! If the CLUBB model is running by itself, but is using a stretched grid
! entered on thermodynamic levels (grid_type = 2), it needs to use the
! thermodynamic level altitudes as input.
! If the CLUBB model is running by itself, but is using a stretched grid
! entered on momentum levels (grid_type = 3), it needs to use the momentum
! level altitudes as input.
real( kind = core_rknd ), intent(in), dimension(gr%nz) :: &
momentum_heights, & ! Momentum level altitudes (input) [m]
thermodynamic_heights ! Thermodynamic level altitudes (input) [m]
integer :: k
! ---- Begin Code ----
if ( .not. l_implemented ) then
if ( grid_type == 1 ) then
! Evenly-spaced grid.
! Momentum level altitudes are defined based on the grid starting
! altitude, zm_init, the constant grid-spacing, deltaz, and the number
! of grid levels, gr%nz.
! Define momentum level altitudes. The first momentum level is at
! altitude zm_init.
do k = 1, gr%nz, 1
gr%zm(k) = zm_init + real( k-1, kind = core_rknd ) * deltaz
enddo
! Define thermodynamic level altitudes. Thermodynamic level altitudes
! are located at the central altitude levels, exactly halfway between
! momentum level altitudes. The lowermost thermodynamic level is
! found by taking 1/2 the altitude difference between the bottom two
! momentum levels and subtracting that value from the bottom momentum
! level. The first thermodynamic level is below zm_init.
gr%zt(1) = zm_init - ( 0.5_core_rknd * deltaz )
do k = 2, gr%nz, 1
gr%zt(k) = 0.5_core_rknd * ( gr%zm(k) + gr%zm(k-1) )
enddo
elseif ( grid_type == 2 ) then
! Stretched (unevenly-spaced) grid: stretched thermodynamic levels.
! Thermodynamic levels are defined according to a stretched grid that
! is entered through the use of an input file. This is similar to a
! SAM-style stretched grid.
! Define thermodynamic level altitudes.
do k = 1, gr%nz, 1
gr%zt(k) = thermodynamic_heights(k)
enddo
! Define momentum level altitudes. Momentum level altitudes are
! located at the central altitude levels, exactly halfway between
! thermodynamic level altitudes. The uppermost momentum level
! altitude is found by taking 1/2 the altitude difference between the
! top two thermodynamic levels and adding that value to the top
! thermodynamic level.
do k = 1, gr%nz-1, 1
gr%zm(k) = 0.5_core_rknd * ( gr%zt(k+1) + gr%zt(k) )
enddo
gr%zm(gr%nz) = gr%zt(gr%nz) + &
0.5_core_rknd * ( gr%zt(gr%nz) - gr%zt(gr%nz-1) )
elseif ( grid_type == 3 ) then
! Stretched (unevenly-spaced) grid: stretched momentum levels.
! Momentum levels are defined according to a stretched grid that is
! entered through the use of an input file. This is similar to a
! WRF-style stretched grid.
! Define momentum level altitudes.
do k = 1, gr%nz, 1
gr%zm(k) = momentum_heights(k)
enddo
! Define thermodynamic level altitudes. Thermodynamic level altitudes
! are located at the central altitude levels, exactly halfway between
! momentum level altitudes. The lowermost thermodynamic level
! altitude is found by taking 1/2 the altitude difference between the
! bottom two momentum levels and subtracting that value from the
! bottom momentum level.
gr%zt(1) = gr%zm(1) - 0.5_core_rknd * ( gr%zm(2) - gr%zm(1) )
do k = 2, gr%nz, 1
gr%zt(k) = 0.5_core_rknd * ( gr%zm(k) + gr%zm(k-1) )
enddo
else
! Invalid grid type.
write(fstderr,*) err_header, "Invalid grid type: ", grid_type, &
". Valid options are 1, 2, or 3."
err_code = clubb_fatal_error
return
endif
else
! The CLUBB parameterization is implemented in a host model.
! Use the host model's momentum level altitudes and thermodynamic level
! altitudes to set up the CLUBB grid.
! Momentum level altitudes from host model.
do k = 1, gr%nz, 1
gr%zm(k) = momentum_heights(k)
enddo
! Thermodynamic level altitudes from host model after possible grid-index
! adjustment for CLUBB interface.
do k = 1, gr%nz, 1
gr%zt(k) = thermodynamic_heights(k)
enddo
endif ! not l_implemented
! Define dzm, the spacing between thermodynamic grid levels; centered over
! momentum grid levels
do k=1,gr%nz-1
gr%dzm(k) = gr%zt(k+1) - gr%zt(k)
enddo
gr%dzm(gr%nz) = gr%dzm(gr%nz-1)
! Define dzt, the spacing between momentum grid levels; centered over
! thermodynamic grid levels
do k=2,gr%nz
gr%dzt(k) = gr%zm(k) - gr%zm(k-1)
enddo
gr%dzt(1) = gr%dzt(2)
! Define invrs_dzm, which is the inverse spacing between thermodynamic grid
! levels; centered over momentum grid levels.
do k=1,gr%nz-1
gr%invrs_dzm(k) = 1._core_rknd / ( gr%zt(k+1) - gr%zt(k) )
enddo
gr%invrs_dzm(gr%nz) = gr%invrs_dzm(gr%nz-1)
! Define invrs_dzt, which is the inverse spacing between momentum grid
! levels; centered over thermodynamic grid levels.
do k=2,gr%nz
gr%invrs_dzt(k) = 1._core_rknd / ( gr%zm(k) - gr%zm(k-1) )
enddo
gr%invrs_dzt(1) = gr%invrs_dzt(2)
! Interpolation Weights: zm grid to zt grid.
! The grid index (k) is the index of the level on the thermodynamic (zt)
! grid. The result is the weights of the upper and lower momentum levels
! (that sandwich the thermodynamic level) applied to that thermodynamic
! level. These weights are normally used in situations where a momentum
! level variable is being solved for implicitly in an equation, and the
! aforementioned variable needs to be interpolated from three successive
! momentum levels (the central momentum level, as well as one momentum level
! above and below the central momentum level) to the intermediate
! thermodynamic grid levels that sandwich the central momentum level.
! For more information, see the comments in function interpolated_aztk_imp.
do k = 1, gr%nz, 1
gr%weights_zm2zt(m_above:m_below,k) &
= interp_weights_zm2zt_imp( k )
enddo
! Interpolation Weights: zt grid to zm grid.
! The grid index (k) is the index of the level on the momentum (zm) grid.
! The result is the weights of the upper and lower thermodynamic levels
! (that sandwich the momentum level) applied to that momentum level. These
! weights are normally used in situations where a thermodynamic level
! variable is being solved for implicitly in an equation, and the
! aforementioned variable needs to be interpolated from three successive
! thermodynamic levels (the central thermodynamic level, as well as one
! thermodynamic level above and below the central thermodynamic level) to
! the intermediate momentum grid levels that sandwich the central
! thermodynamic level.
! For more information, see the comments in function interpolated_azmk_imp.
do k = 1, gr%nz, 1
gr%weights_zt2zm(t_above:t_below,k) &
= interp_weights_zt2zm_imp( k )
enddo
return
end subroutine setup_grid_heights
!=============================================================================
subroutine read_grid_heights( nzmax, grid_type, &
zm_grid_fname, zt_grid_fname, &
file_unit, &
momentum_heights, &
thermodynamic_heights )
! Description:
! This subroutine is used foremost in cases where the grid_type corresponds
! with the stretched (unevenly-spaced) grid options (either grid_type = 2 or
! grid_type = 3). This subroutine reads in the values of the stretched grid
! altitude levels for either the thermodynamic level grid or the momentum
! level grid. This subroutine also handles basic error checking for all
! three grid types.
!------------------------------------------------------------------------
use constants_clubb, only: &
fstderr ! Variable(s)
use file_functions, only: &
file_read_1d ! Procedure(s)
use clubb_precision, only: &
core_rknd ! Variable(s)
use error_code, only: &
err_code, & ! Error indicator
clubb_fatal_error, & ! Constant
err_header ! String
implicit none
! Input Variables.
! Declared number of vertical levels.
integer, intent(in) :: &
nzmax
! If CLUBB is running on it's own, this option determines if it is using:
! 1) an evenly-spaced grid;
! 2) a stretched (unevenly-spaced) grid entered on the thermodynamic grid
! levels (with momentum levels set halfway between thermodynamic levels);
! or
! 3) a stretched (unevenly-spaced) grid entered on the momentum grid levels
! (with thermodynamic levels set halfway between momentum levels).
integer, intent(in) :: &
grid_type
character(len=*), intent(in) :: &
zm_grid_fname, & ! Path and filename of file for momentum level altitudes
zt_grid_fname ! Path and filename of file for thermodynamic level altitudes
integer, intent(in) :: &
file_unit ! Unit number for zt_grid_fname & zm_grid_fname (based on the OpenMP thread)
! Output Variables.
! If the CLUBB model is running by itself, but is using a stretched grid
! entered on thermodynamic levels (grid_type = 2), it needs to use the
! thermodynamic level altitudes as input.
! If the CLUBB model is running by itself, but is using a stretched grid
! entered on momentum levels (grid_type = 3), it needs to use the momentum
! level altitudes as input.
real( kind = core_rknd ), dimension(nzmax), intent(out) :: &
momentum_heights, & ! Momentum level altitudes (file input) [m]
thermodynamic_heights ! Thermodynamic level altitudes (file input) [m]
! Local Variables.
integer :: &
zt_level_count, & ! Number of altitudes found in zt_grid_fname
zm_level_count ! Number of altitudes found in zm_grid_fname
integer :: input_status ! Status of file being read:
! > 0 ==> error reading file.
! = 0 ==> no error and more file to be read.
! < 0 ==> end of file indicator.
! Generic variable for storing file data while counting the number
! of file entries.
real( kind = core_rknd ) :: generic_input_item
integer :: k ! Loop index
! ---- Begin Code ----
! Declare the momentum level altitude array and the thermodynamic level
! altitude array to be 0 until overwritten.
momentum_heights(1:nzmax) = 0.0_core_rknd
thermodynamic_heights(1:nzmax) = 0.0_core_rknd
! Avoid uninitialized memory
generic_input_item = 0.0_core_rknd
if ( grid_type == 1 ) then
! Evenly-spaced grid.
! Grid level altitudes are based on a constant distance between them and
! a starting point for the bottom of the grid.
! As a way of error checking, make sure that there isn't any file entry
! for either momentum level altitudes or thermodynamic level altitudes.
if ( zm_grid_fname /= '' ) then
write(fstderr,*) err_header, &
"An evenly-spaced grid has been selected. " &
// " Please reset zm_grid_fname to ''."
err_code = clubb_fatal_error
return
endif
if ( zt_grid_fname /= '' ) then
write(fstderr,*) err_header, &
"An evenly-spaced grid has been selected. " &
// " Please reset zt_grid_fname to ''."
err_code = clubb_fatal_error
return
endif
elseif ( grid_type == 2 ) then
! Stretched (unevenly-spaced) grid: stretched thermodynamic levels.
! Thermodynamic levels are defined according to a stretched grid that is
! entered through the use of an input file. Momentum levels are set
! halfway between thermodynamic levels. This is similar to a SAM-style
! stretched grid.
! As a way of error checking, make sure that there isn't any file entry
! for momentum level altitudes.
if ( zm_grid_fname /= '' ) then
write(fstderr,*) err_header, &
"Thermodynamic level altitudes have been selected " &
// "for use in a stretched (unevenly-spaced) grid. " &
// " Please reset zm_grid_fname to ''."
err_code = clubb_fatal_error
return
endif
!$omp critical
! Open the file zt_grid_fname.
open( unit=file_unit, file=zt_grid_fname, &
status='old', action='read' )
! Find the number of thermodynamic level altitudes listed
! in file zt_grid_fname.
zt_level_count = 0
do
read( unit=file_unit, fmt=*, iostat=input_status ) &
generic_input_item
if ( input_status < 0 ) exit ! end of file indicator
if ( input_status > 0 ) then
write(fstderr,*) err_header, & ! error reading input
"Error reading thermodynamic level input file."
err_code = clubb_fatal_error
exit
end if
zt_level_count = zt_level_count + 1
enddo
! Close the file zt_grid_fname.
close( unit=file_unit )
!$omp end critical
if ( err_code == clubb_fatal_error ) return
! Check that the number of thermodynamic grid altitudes in the input file
! matches the declared number of CLUBB grid levels (nzmax).
if ( zt_level_count /= nzmax ) then
write(fstderr,*) &
"The number of thermodynamic grid altitudes " &
// "listed in file " // trim(zt_grid_fname) &
// " does not match the number of CLUBB grid " &
// "levels specified in the model.in file."
write(fstderr,*) &
"Number of thermodynamic grid altitudes listed: ", &
zt_level_count
write(fstderr,*) &
"Number of CLUBB grid levels specified: ", nzmax
err_code = clubb_fatal_error
return
endif
! Read the thermodynamic level altitudes from zt_grid_fname.
call file_read_1d( file_unit, zt_grid_fname, nzmax, 1, &
thermodynamic_heights )
! Check that each thermodynamic level altitude increases
! in height as the thermodynamic level grid index increases.
do k = 2, nzmax, 1
if ( thermodynamic_heights(k) &
<= thermodynamic_heights(k-1) ) then
write(fstderr,*) &
"The declared thermodynamic level grid " &
// "altitudes are not increasing in height " &
// "as grid level index increases."
write(fstderr,*) &
"Grid index: ", k-1, ";", &
" Thermodynamic level altitude: ", &
thermodynamic_heights(k-1)
write(fstderr,*) &
"Grid index: ", k, ";", &
" Thermodynamic level altitude: ", &
thermodynamic_heights(k)
err_code = clubb_fatal_error
return
endif
enddo
elseif ( grid_type == 3 ) then
! Stretched (unevenly-spaced) grid: stretched momentum levels.