-
Notifications
You must be signed in to change notification settings - Fork 11
/
LDAS_Forcing.F90
5572 lines (3888 loc) · 192 KB
/
LDAS_Forcing.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
#include "MAPL_Generic.h"
module LDAS_ForceMod
! collection of *forcing* subroutines for enkf_driver
! (originally these routines were in clsm_ensdrv_drv_routines.F90)
! reichle, 13 Aug 2008
use, intrinsic :: iso_c_binding
use ESMF
use MAPL_Mod
use MAPL_ShmemMod
use LDAS_ensdrv_Globals, ONLY: &
logunit, &
logit, &
root_logit, &
nodata_generic, &
nodata_tol_generic, &
nodata_tolfrac_generic
use MAPL_ConstantsMod, ONLY: &
stefan_boltzmann => MAPL_STFBOL, &
Tzero => MAPL_TICE
use MAPL_SatVaporMod, ONLY: &
MAPL_EQsat
use LDAS_DriverTypes, ONLY: &
met_force_type
use LDAS_TileCoordType, ONLY: &
tile_coord_type
use LDAS_DateTimeMod, ONLY: &
date_time_type, &
augment_date_time, &
datetime_lt_refdatetime, &
datetime_le_refdatetime, &
is_leap_year, &
get_dofyr_pentad
use LDAS_ExceptionsMod, ONLY: &
ldas_abort, &
LDAS_GENERIC_ERROR
use RepairForcingMod, ONLY: &
repair_forcing
use LDAS_HashTable, only: &
Hash_Table
implicit none
include 'netcdf.inc'
! everything is private by default unless made public
private
public :: get_forcing
public :: LDAS_move_new_force_to_old
public :: GEOS_closefile
type(Hash_Table),public :: FileOpenedHash
!public :: ignore_SWNET_for_snow
!logical :: ignore_SWNET_for_snow ! fixes sibalb bug in MERRA snow albedo
real, parameter :: DEFAULT_REFH = 10. ! m
! real, parameter :: SWDN_MAX = 1360. ! W/m2
! real, parameter :: LWDN_EMISS_MIN = 0.5 ! min effective emissivity for LWdown
! real, parameter :: LWDN_EMISS_MAX = 1.0 ! max effective emissivity for LWdown
! The above three parameters are moved to repairForcingmod ! W.Jiang
character(10), private :: tmpstring10
character(40), private :: tmpstring40
real, contiguous, pointer :: ptrShForce(:,:)=>null()
type local_grid
integer :: N_lon = 0
integer :: N_lat = 0
integer :: N_cat = 0
integer,allocatable :: i1(:),i2(:),j1(:),j2(:)
real,allocatable :: x1(:),x2(:),y1(:),y2(:)
end type local_grid
type(local_grid), target :: local_info
contains
! ********************************************************************
subroutine get_forcing( date_time, force_dtstep, met_path, met_tag, &
N_catd, tile_coord, met_hinterp, &
MERRA_file_specs, GEOS_Forcing, met_force_obs_tile_new, &
AEROSOL_DEPOSITION,init, alb_from_SWnet )
! Read and check meteorological forcing data for the domain.
!
! time convention:
! - forcing states (such as Tair) are snapshots at date_time
! - forcing fluxes (such as SWdn) are time avg over *subsequent* forcing
! interval (date_time:date_time+force_dtstep)
!
! The above time convention is heritage from older versions of the
! off-line driver and creates problems with "operational" forcing
! data from GEOS5. For "operational" integrations, the forward-looking
! forcing fluxes are not available for "met_force_obs_tile_new".
!
! As a work-around, the output parameter "move_met_force_obs_new_to_old"
! is used to treat "operational" forcing data sets accordingly in the
! main program. This work-around replaces an older work-around that
! was less efficient.
!
! When LDASsa is integrated within the coupled GEOS5 DAS, initial (time-avg)
! "tavg1_2d_*_Nx" files are not available. Use optional "init" flag to
! deal with this situation.
!
! reichle, 28 March 2006
! reichle, 13 March 2008 - added optional "init" flag
! qliu+reichle, 12 Aug 2008 - new field RefH (reference height) in met_force_type
! reichle, 25 Sep 2009 - removed unneeded inputs
! reichle, 23 Feb 2016 - new and more efficient work-around to make GEOS-5
! forcing work with LDASsa time convention for forcing data
implicit none
type(date_time_type), intent(in) :: date_time
integer, intent(in) :: force_dtstep
character(*), intent(in) :: met_path
character(*), intent(in) :: met_tag
integer, intent(in) :: N_catd, met_hinterp
type(tile_coord_type), dimension(:), pointer :: tile_coord ! input
logical,intent(out) :: MERRA_file_specs
logical,intent(out) :: GEOS_forcing
type(met_force_type), dimension(N_catd), intent(out) :: &
met_force_obs_tile_new
integer,intent(in) :: AEROSOL_DEPOSITION
logical, intent(in), optional :: init, alb_from_SWnet
! local variables
real :: nodata_forcing
type(date_time_type) :: date_time_tmp
character(len=*), parameter :: Iam = 'get_forcing'
character(len=400) :: err_msg
! --------------------------------------------------------------
!
! shift forcing date if so indicated by met_tag (for twin experiments,
! see function shift_forcing_date for details) - reichle, 6 Apr 2007
date_time_tmp = shift_forcing_date(met_tag, date_time)
! set reference height to default value (if appropriate, will be overwriten
! within specific subroutine)
met_force_obs_tile_new%RefH = DEFAULT_REFH
! set SWnet, PARdrct, PARdffs to nodata_generic
! (Note that nodata_forcing is set to the native no-data-value
! in the individual get_*() subroutines and used to communicate with
! check_forcing_nodata. AFTER the call to check_forcing_nodata all forcing
! fields EXCEPT SWnet must NOT be no-data values, and SWnet must be
! nodata_generic if unavailable.)
!
! reichle+qliu, 8 Oct 2008
! reichle, 23 Feb 2009 -- same goes for ParDrct, ParDffs
! reichle, 5 Mar 2009 -- deleted ParDrct, ParDffs after testing found no impact
! reichle, 22 Jul 2010 -- fixed treatment of SWnet nodata values
! reichle, 20 Dec 2011 -- reinstated PARdrct and PARdffs for MERRA-Land file specs
! ---------------------------------------------------------------------------------
!
! initialize
MERRA_file_specs = .false.
GEOS_forcing = .false.
met_force_obs_tile_new%SWnet = nodata_generic
met_force_obs_tile_new%PARdrct = nodata_generic
met_force_obs_tile_new%PARdffs = nodata_generic
! ---------------------------------------------------------------------------------
!
! get forcing in tile space
if (index(met_tag, 'Berg_netcdf')/=0) then
call get_Berg_netcdf( date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
elseif (index(met_tag, 'GLDAS_2x2_5_netcdf')/=0) then
call get_GLDAS_2x2_5_netcdf(date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
elseif (index(met_tag, 'Viviana_OK')/=0) then
! vmaggion & reichle, 17 July 2008
!
! use 2x2.5 deg GLDAS for all forcing fields except precip
call get_GLDAS_2x2_5_netcdf(date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
if (index(met_tag, 'Viviana_OK_nopert')/=0) then
call get_Viviana_OK_precip(10, date_time_tmp, met_path, met_tag, &
N_catd, tile_coord, met_force_obs_tile_new)
end if
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
elseif (index(met_tag, 'GSWP2_1x1_netcdf')/=0) then
call get_GSWP2_1x1_netcdf( date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
elseif (index(met_tag, 'RedArk_ASCII')/=0) then
call get_RedArk_ASCII( date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
elseif (index(met_tag, 'RedArk_GOLD')/=0) then
call get_RedArk_GOLD( date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
elseif (index(met_tag, 'RedArk_Princeton')/=0) then
call get_RedArk_Princeton( date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
elseif (index(met_tag, 'Princeton_netcdf')/=0) then ! tyamada+reichle, 17 Jul 2007
call get_Princeton_netcdf( date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
elseif (index(met_tag, 'conus_0.5d_netcdf')/=0) then ! sarith+reichle, 17 Jul 2007
call get_conus_netcdf( date_time_tmp, met_path, N_catd, tile_coord, &
met_force_obs_tile_new, nodata_forcing)
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord )
else ! assume forcing from GEOS5 GCM ("DAS" or "MERRA") output
if(root_logit) write (logunit,*) 'get_forcing(): assuming GEOS-5 forcing data set'
GEOS_forcing = .true.
! note "met_tag" in call to get_GEOSgcm_gfio (interface differs
! from other get_* subroutines)
!call get_GEOSgcm_gfio( date_time_tmp, met_path, met_tag, &
! N_catd, tile_coord, &
! met_force_obs_tile_new, nodata_forcing)
call get_GEOS( date_time_tmp, force_dtstep, &
met_path, met_tag, N_catd, tile_coord, met_hinterp, &
met_force_obs_tile_new, nodata_forcing, MERRA_file_specs, &
AEROSOL_DEPOSITION, init )
! check for nodata values and unphysical values
! (check here, not outside "if" block, because of GEOSgcm case)
call check_forcing_nodata( N_catd, tile_coord, nodata_forcing, &
met_force_obs_tile_new )
! call repair_forcing with switch "unlimited_Qair=.true." for GEOS5 forcing
! (default is to limit Qair so that it does not exceed Qair_sat)
! reichle+qliu, 8 Oct 2008
!
! likewise for "unlimited_LWdown=.true."
! reichle, 11 Feb 2009
call repair_forcing( N_catd, met_force_obs_tile_new, &
echo=.true., tile_coord=tile_coord, &
unlimited_Qair=.true., unlimited_LWdown=.true. )
! Subroutine get_GEOS() reads forcing fluxes from "previous"
! interval, not from "subsequent" interval, because in operational
! applications the "subsequent" fluxes for "met_force_new" are not
! available. The following lines restore consistency with the
! time convention stated above. Note that only "old" fluxes
! are needed in subroutine interpolate_to_timestep(), and
! "met_force_obs_tile_new" is set to nodata for forcing fluxes.
! The calls below are moved to LDAS_move_new_force_to_old
! met_force_obs_tile_old%Rainf_C = met_force_obs_tile_new%Rainf_C
! met_force_obs_tile_old%Rainf = met_force_obs_tile_new%Rainf
! met_force_obs_tile_old%Snowf = met_force_obs_tile_new%Snowf
! met_force_obs_tile_old%LWdown = met_force_obs_tile_new%LWdown
! met_force_obs_tile_old%SWdown = met_force_obs_tile_new%SWdown
! met_force_obs_tile_old%SWnet = met_force_obs_tile_new%SWnet
! met_force_obs_tile_old%PARdrct = met_force_obs_tile_new%PARdrct
! met_force_obs_tile_old%PARdffs = met_force_obs_tile_new%PARdffs
! treat Wind as flux when forcing with MERRA
! if (MERRA_file_specs) met_force_obs_tile_old%Wind = met_force_obs_tile_new%Wind
! met_force_obs_tile_new%Rainf_C = nodata_generic
! met_force_obs_tile_new%Rainf = nodata_generic
! met_force_obs_tile_new%Snowf = nodata_generic
! met_force_obs_tile_new%LWdown = nodata_generic
! met_force_obs_tile_new%SWdown = nodata_generic
! met_force_obs_tile_new%SWnet = nodata_generic
! met_force_obs_tile_new%PARdrct = nodata_generic
! met_force_obs_tile_new%PARdffs = nodata_generic
!
! if (MERRA_file_specs) met_force_obs_tile_new%Wind = nodata_generic
end if
! make sure SWnet is generally available if needed to back out albedo
! (only works for GEOS forcing, e.g., MERRA, FP, FP-IT)
!
! NOTE: need to check here because GEOS forcing is the default
! forcing data set in the "if... elseif... elseif... else..."
! statement above, that is, it is only asserted here whether
! the requested forcing data are GEOS.
if (present(alb_from_SWnet)) then
if ( (.not. GEOS_forcing) .and. (alb_from_SWnet) ) then
! stop if per nml inputs the albedo should be backed
! out from SWnet but forcing data is not GEOS
err_msg = 'requested nml input [alb_from_SWnet=.true.] ' // &
'*only* works with GEOS forcing'
call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg)
end if
end if
! stop if anything other than nearest-neighbor met forcing interpolation was
! requested for a non-GEOS5 forcing dataset
if ((.not. GEOS_forcing) .and. (met_hinterp>0)) then
err_msg = 'for non-GEOS forcing, only ' // &
'nearest-neighbor interpolation is available'
call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg)
end if
end subroutine get_forcing
!*****************************
subroutine LDAS_move_new_force_to_old(new_force,old_force,MERRA_file_specs, GEOS_forcing,AEROSOL_DEPOSITION)
implicit none
type(met_force_type), dimension(:), intent(inout) :: new_force
type(met_force_type), dimension(:), intent(inout) :: old_force
logical,intent(in) :: MERRA_file_specs
logical,intent(in) :: GEOS_forcing
integer, intent(in) :: AEROSOL_DEPOSITION
old_force%Rainf_C = 0.0
old_force%Rainf = 0.0
old_force%Snowf = 0.0
old_force%LWdown = 0.0
old_force%SWdown = 0.0
old_force%SWnet = 0.0
old_force%PARdrct = 0.0
old_force%PARdffs = 0.0
if (.not. GEOS_forcing) return
old_force%Rainf_C = new_force%Rainf_C
old_force%Rainf = new_force%Rainf
old_force%Snowf = new_force%Snowf
old_force%LWdown = new_force%LWdown
old_force%SWdown = new_force%SWdown
old_force%SWnet = new_force%SWnet
old_force%PARdrct = new_force%PARdrct
old_force%PARdffs = new_force%PARdffs
new_force%Rainf_C = nodata_generic
new_force%Rainf = nodata_generic
new_force%Snowf = nodata_generic
new_force%LWdown = nodata_generic
new_force%SWdown = nodata_generic
new_force%SWnet = nodata_generic
new_force%PARdrct = nodata_generic
new_force%PARdffs = nodata_generic
if( AEROSOL_DEPOSITION /=0) then
old_force%DUDP001 = new_force%DUDP001
old_force%DUDP002 = new_force%DUDP002
old_force%DUDP003 = new_force%DUDP003
old_force%DUDP004 = new_force%DUDP004
old_force%DUDP005 = new_force%DUDP005
old_force%DUSV001 = new_force%DUSV001
old_force%DUSV002 = new_force%DUSV002
old_force%DUSV003 = new_force%DUSV003
old_force%DUSV004 = new_force%DUSV004
old_force%DUSV005 = new_force%DUSV005
old_force%DUWT001 = new_force%DUWT001
old_force%DUWT002 = new_force%DUWT002
old_force%DUWT003 = new_force%DUWT003
old_force%DUWT004 = new_force%DUWT004
old_force%DUWT005 = new_force%DUWT005
old_force%DUSD001 = new_force%DUSD001
old_force%DUSD002 = new_force%DUSD002
old_force%DUSD003 = new_force%DUSD003
old_force%DUSD004 = new_force%DUSD004
old_force%DUSD005 = new_force%DUSD005
old_force%BCDP001 = new_force%BCDP001
old_force%BCDP002 = new_force%BCDP002
old_force%BCSV001 = new_force%BCSV001
old_force%BCSV002 = new_force%BCSV002
old_force%BCWT001 = new_force%BCWT001
old_force%BCWT002 = new_force%BCWT002
old_force%BCSD001 = new_force%BCSD001
old_force%BCSD002 = new_force%BCSD002
old_force%OCDP001 = new_force%OCDP001
old_force%OCDP002 = new_force%OCDP002
old_force%OCSV001 = new_force%OCSV001
old_force%OCSV002 = new_force%OCSV002
old_force%OCWT001 = new_force%OCWT001
old_force%OCWT002 = new_force%OCWT002
old_force%OCSD001 = new_force%OCSD001
old_force%OCSD002 = new_force%OCSD002
old_force%SUDP003 = new_force%SUDP003
old_force%SUSV003 = new_force%SUSV003
old_force%SUWT003 = new_force%SUWT003
old_force%SUSD003 = new_force%SUSD003
old_force%SSDP001 = new_force%SSDP001
old_force%SSDP002 = new_force%SSDP002
old_force%SSDP003 = new_force%SSDP003
old_force%SSDP004 = new_force%SSDP004
old_force%SSDP005 = new_force%SSDP005
old_force%SSSV001 = new_force%SSSV001
old_force%SSSV002 = new_force%SSSV002
old_force%SSSV003 = new_force%SSSV003
old_force%SSSV004 = new_force%SSSV004
old_force%SSSV005 = new_force%SSSV005
old_force%SSWT001 = new_force%SSWT001
old_force%SSWT002 = new_force%SSWT002
old_force%SSWT003 = new_force%SSWT003
old_force%SSWT004 = new_force%SSWT004
old_force%SSWT005 = new_force%SSWT005
old_force%SSSD001 = new_force%SSSD001
old_force%SSSD002 = new_force%SSSD002
old_force%SSSD003 = new_force%SSSD003
old_force%SSSD004 = new_force%SSSD004
old_force%SSSD005 = new_force%SSSD005
new_force%DUDP001 = nodata_generic
new_force%DUDP002 = nodata_generic
new_force%DUDP003 = nodata_generic
new_force%DUDP004 = nodata_generic
new_force%DUDP005 = nodata_generic
new_force%DUSV001 = nodata_generic
new_force%DUSV002 = nodata_generic
new_force%DUSV003 = nodata_generic
new_force%DUSV004 = nodata_generic
new_force%DUSV005 = nodata_generic
new_force%DUWT001 = nodata_generic
new_force%DUWT002 = nodata_generic
new_force%DUWT003 = nodata_generic
new_force%DUWT004 = nodata_generic
new_force%DUWT005 = nodata_generic
new_force%DUSD001 = nodata_generic
new_force%DUSD002 = nodata_generic
new_force%DUSD003 = nodata_generic
new_force%DUSD004 = nodata_generic
new_force%DUSD005 = nodata_generic
new_force%BCDP001 = nodata_generic
new_force%BCDP002 = nodata_generic
new_force%BCSV001 = nodata_generic
new_force%BCSV002 = nodata_generic
new_force%BCWT001 = nodata_generic
new_force%BCWT002 = nodata_generic
new_force%BCSD001 = nodata_generic
new_force%BCSD002 = nodata_generic
new_force%OCDP001 = nodata_generic
new_force%OCDP002 = nodata_generic
new_force%OCSV001 = nodata_generic
new_force%OCSV002 = nodata_generic
new_force%OCWT001 = nodata_generic
new_force%OCWT002 = nodata_generic
new_force%OCSD001 = nodata_generic
new_force%OCSD002 = nodata_generic
new_force%SUDP003 = nodata_generic
new_force%SUSV003 = nodata_generic
new_force%SUWT003 = nodata_generic
new_force%SUSD003 = nodata_generic
new_force%SSDP001 = nodata_generic
new_force%SSDP002 = nodata_generic
new_force%SSDP003 = nodata_generic
new_force%SSDP004 = nodata_generic
new_force%SSDP005 = nodata_generic
new_force%SSSV001 = nodata_generic
new_force%SSSV002 = nodata_generic
new_force%SSSV003 = nodata_generic
new_force%SSSV004 = nodata_generic
new_force%SSSV005 = nodata_generic
new_force%SSWT001 = nodata_generic
new_force%SSWT002 = nodata_generic
new_force%SSWT003 = nodata_generic
new_force%SSWT004 = nodata_generic
new_force%SSWT005 = nodata_generic
new_force%SSSD001 = nodata_generic
new_force%SSSD002 = nodata_generic
new_force%SSSD003 = nodata_generic
new_force%SSSD004 = nodata_generic
new_force%SSSD005 = nodata_generic
endif ! AEROSOL_DEPOSITION /=0
! treat Wind as flux when forcing with MERRA
if (MERRA_file_specs) then
old_force%Wind = new_force%Wind
new_force%Wind = nodata_generic
endif
end subroutine LDAS_move_new_force_to_old
! ****************************************************************
subroutine get_Berg_netcdf(date_time, met_path, N_catd, tile_coord, &
met_force_new, nodata_forcing )
! read Berg_NetCDF files and extract forcings in tile space
! (uses nearest neighbor interpolation)
! reichle, 25 May 2005
! reichle, 23 Feb 2009 -- revised treatment of Rainf_C
implicit none
type(date_time_type), intent(in) :: date_time
character(*), intent(in) :: met_path
integer, intent(in) :: N_catd
type(tile_coord_type), dimension(:), pointer :: tile_coord ! input
type(met_force_type), dimension(N_catd), intent(inout) :: met_force_new
real, intent(out) :: nodata_forcing
! Berg grid and netcdf parameters
integer, parameter :: berg_grid_N_lon = 720
integer, parameter :: berg_grid_N_lat = 360
real, parameter :: berg_grid_ll_lon = -180.
real, parameter :: berg_grid_ll_lat = -90.
real, parameter :: berg_grid_dlon = .5
real, parameter :: berg_grid_dlat = .5
integer, parameter :: N_berg_compressed = 67420
! Berg forcing time step in hours
integer, parameter :: dt_berg_in_hours = 6
integer, parameter :: nciv_land_i = 3
integer, parameter :: nciv_land_j = 4
integer, parameter :: nciv_data = 7
integer, parameter :: N_berg_vars = 8
real, parameter :: nodata_berg = 1.e20
character(40), dimension(N_berg_vars), parameter :: berg_dir = (/ &
'PREC-CONV/', &
'PREC-TOTL/', &
'PRES-SRF/ ', &
'RAD-LW/ ', &
'RAD-SW/ ', &
'TEMP-AIR/ ', &
'TEMP-DEW/ ', &
'WIND/ ' /)
character(40), dimension(N_berg_vars), parameter :: berg_name = (/ &
'RainfSnowf_C_ecmwf', &
'RainfSnowf_ecmwf ', &
'PSurf_ecmwf ', &
'LWdown_ecmwf ', &
'SWdown_ecmwf ', &
'Tair_ecmwf ', &
'Tdew_ecmwf ', &
'Wind_ecmwf ' /)
! local variables
!!real, parameter :: Tzero = 273.16
real :: tol
real, dimension(berg_grid_N_lon,berg_grid_N_lat) :: tmp_grid
integer, dimension(N_berg_compressed) :: land_i_berg, land_j_berg
integer, dimension(N_catd) :: i_ind, j_ind
real, dimension(N_berg_compressed) :: tmp_vec
real, dimension(N_catd,N_berg_vars) :: force_array
integer, dimension(2) :: start, icount
integer :: k, n, hours_in_month, berg_var, ierr, ncid
real :: this_lon, this_lat, dt_berg_in_seconds
character(4) :: YYYY
character(2) :: MM
character(300) :: fname
character(len=*), parameter :: Iam = 'get_Berg_netcdf'
character(len=400) :: err_msg
! --------------------------------------------------------------------
dt_berg_in_seconds = real(3600*dt_berg_in_hours)
nodata_forcing = nodata_berg
tol = abs(nodata_forcing*nodata_tolfrac_generic)
! assemble year and month strings
write (YYYY,'(i4.4)') date_time%year
write (MM, '(i2.2)') date_time%month
! find out which data are needed
! compressed space dimension (always read global vector)
start(1) = 1
icount(1) = N_berg_compressed
! time dimension (first entry in Berg_NetCDF file is at 0Z)
if ( (date_time%min/=0) .or. (date_time%sec/=0) .or. &
(mod(date_time%hour,dt_berg_in_hours)/=0) ) then
call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'timing error')
end if
hours_in_month = (date_time%day-1)*24 + date_time%hour
start(2) = hours_in_month / dt_berg_in_hours + 1
icount(2) = 1
! ----------------------------------------------
!
! compute indices for nearest neighbor interpolation from Berg grid
! to tile space
!
! (NOTE: this should at some point be replaced with a regridding
! subroutine that interpolates from the
! native forcing grid to the GCM atmospheric grid that is used
! to cut catchments into tiles - then "standard" grid2tile
! using tile_coord%atm_i and tile_coord%atm_j applies.
! reichle, 26 May 2005)
do k=1,N_catd
! ll_lon and ll_lat refer to lower left corner of grid cell
! (as opposed to the grid point in the center of the grid cell)
this_lon = tile_coord(k)%com_lon
this_lat = tile_coord(k)%com_lat
i_ind(k) = ceiling( (this_lon - berg_grid_ll_lon)/berg_grid_dlon )
j_ind(k) = ceiling( (this_lat - berg_grid_ll_lat)/berg_grid_dlat )
! NOTE: For a "date line on center" grid and (180-dlon/2) < lon < 180
! we now have i_ind=(grid%N_lon+1)
! This would need to be fixed.
if (i_ind(k)>berg_grid_N_lon) i_ind(k)=1
end do
! ------------------------------------------------------
!
! read compression parameters (same for all data variables and time steps)
berg_var = 1
fname = trim(met_path) // trim(berg_dir(berg_var)) // '/' // YYYY &
// '/' // trim(berg_name(berg_var)) // '.' // YYYY // MM // '.nc'
if(root_logit) write(logunit,*) 'get netcdf compression params from ' // trim(fname)
ierr = NF_OPEN(fname,NF_NOWRITE,ncid)
if (ierr/=0) then
err_msg = 'error opening netcdf file'
call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg)
end if
ierr = NF_GET_VARA_INT(ncid, nciv_land_i, start, icount, land_i_berg)
ierr = NF_GET_VARA_INT(ncid, nciv_land_j, start, icount, land_j_berg)
ierr = NF_CLOSE(ncid)
! ------------------------------------------------------
!
! get forcing data
do berg_var = 1,N_berg_vars
! open file, read compressed data, and put on global grid
fname = trim(met_path) // trim(berg_dir(berg_var)) // '/' // YYYY &
// '/' // trim(berg_name(berg_var)) // '.' // YYYY // MM // '.nc'
if(root_logit) write (logunit,*) 'opening ' // trim(fname)
ierr = NF_OPEN(fname,NF_NOWRITE,ncid)
if (ierr/=0) then
err_msg = 'error opening netcdf file'
call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg)
end if
ierr = NF_GET_VARA_REAL(ncid, nciv_data, start, icount, tmp_vec )
ierr = NF_CLOSE(ncid)
tmp_grid = nodata_forcing
do n=1,N_berg_compressed
tmp_grid(land_i_berg(n), land_j_berg(n) ) = tmp_vec(n)
end do
! interpolate to tile space
! (NOTE: This should at some point be replaced with a regridding
! subroutine that interpolates from the
! native forcing grid to the GCM atmospheric grid that is used
! to cut catchments into tiles - then "standard" grid2tile
! using tile_coord%atm_i and tile_coord%atm_j applies.
! reichle, 26 May 2005)
do k=1,N_catd
force_array(k,berg_var) = tmp_grid(i_ind(k), j_ind(k))
end do
end do
! convert variables and units of force_array to match met_force_type,
! put into structure
! from Berg files:
!
! force_array(:,1) = PREC-CONV = Rainf_C+Snowf_C kg/m2 (6h total) ???
! force_array(:,2) = PREC-TOTL = RainfSnowf kg/m2 (6h total)
! force_array(:,3) = PRES-SRF = PSurf Pa
! force_array(:,4) = RAD-LW = LWdown W/m2
! force_array(:,5) = RAD-SW = SWdown W/m2
! force_array(:,6) = TEMP-AIR = Tair K
! force_array(:,7) = TEMP-DEW = Tdew K
! force_array(:,8) = WIND = Wind m/s
met_force_new%Psurf = force_array(:,3)
met_force_new%LWdown = force_array(:,4)
met_force_new%SWdown = force_array(:,5)
met_force_new%Tair = force_array(:,6)
met_force_new%Wind = force_array(:,8)
! get specific humidity from dew point temperature
do k=1,N_catd
if ( abs(force_array(k,3)-nodata_berg)<tol .or. &
abs(force_array(k,7)-nodata_berg)<tol ) then
met_force_new(k)%Qair = nodata_forcing
else
met_force_new(k)%Qair = MAPL_EQsat(force_array(k,7),force_array(k,3))
end if
end do
! rain and snow:
! convert from Berg 6h totals [kg/m2] (or [mm]) into precipitation
! rates [kg/m2/s]
met_force_new%Rainf = 0.
met_force_new%Rainf_C = 0.
met_force_new%Snowf = 0.
do k=1,N_catd
! set convective precip to zero for no-data-values
if (abs(force_array(k,1)-nodata_berg)<tol) force_array(k,1) = 0.
if (abs(force_array(k,2)-nodata_berg)<tol) then
met_force_new(k)%Snowf = nodata_forcing
met_force_new(k)%Rainf = nodata_forcing
met_force_new(k)%Rainf_C = nodata_forcing
else
if ( met_force_new(k)%Tair < Tzero ) then
met_force_new(k)%Snowf = force_array(k,2)/dt_berg_in_seconds
else
met_force_new(k)%Rainf = force_array(k,2)/dt_berg_in_seconds
! leave Rainf_C=0 until sure that "PREC-CONV = Rainf_C+Snowf_C"
!!met_force_new(k)%Rainf_C = force_array(k,1)/dt_berg_in_seconds
end if
end if
end do
end subroutine get_Berg_netcdf
! ****************************************************************
subroutine get_RedArk_ASCII(date_time, met_path, N_catd, tile_coord, &
met_force_new, nodata_forcing )
! read RedArk ASCII forcing files from Hatim Sharif and Wade Crow
! extract forcings in tile space (use nearest neighbor interpolation)
! reichle, 13 Apr 2006
implicit none
type(date_time_type), intent(in) :: date_time
character(*), intent(in) :: met_path
integer, intent(in) :: N_catd
type(tile_coord_type), dimension(:), pointer :: tile_coord ! input
type(met_force_type), dimension(N_catd), intent(inout) :: met_force_new
real, intent(out) :: nodata_forcing
! RedArk parameters
integer, parameter :: N_redark_basins = 314
character(40), parameter :: fname_b2t = 'RedArkBasins_2_FV288x181Tiles.dat'
! RedArk forcing time step in hours
integer, parameter :: dt_redark_in_hours = 1
integer, parameter :: N_redark_vars = 6
! not checking for "bad" forcing data from Hatim/Wade, but need
! to provide output "nodata_forcing"
real, parameter :: nodata_redark = 1.e20
real, parameter :: RedArk_Psurf_in_Pa = 101300.
! local variables
type(date_time_type) :: date_time_CST
!!real, parameter :: Tzero = 273.16
integer, dimension(N_catd) :: b_ind
real, dimension(N_redark_basins,N_redark_vars) :: force_array
integer :: i, j, k, tmp_tile_id
real :: dt_redark_in_seconds, Tdew
character(4) :: YYYY
character(3) :: DDD
character(2) :: HH
character(300) :: fname
character(len=*), parameter :: Iam = 'get_RedArk_ASCII'
character(len=400) :: err_msg
! --------------------------------------------------------------------
nodata_forcing = nodata_redark
dt_redark_in_seconds = real(3600*dt_redark_in_hours)
! convert GMT (date_time) into CST (Central Standard Time)
date_time_CST = date_time
call augment_date_time( -6*3600, date_time_CST )
! assemble year, day-of-year, and hour strings