-
Notifications
You must be signed in to change notification settings - Fork 0
/
Nest_mod.f90
executable file
·1668 lines (1518 loc) · 68 KB
/
Nest_mod.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 Nest_mod
! This module performs the reading or writing of data for nested runs
!
! The Nesting modes and related setings read from Nest_config nml:
! MODE_READ
! 'NONE': do not read (default)
! 'NHOUR': read every NHOURREAD
! 'START': read at the start of run
! 'RESTART: read at the start of run, if the files are found
! MODE_SAVE
! 'NONE': do not write (default)
! 'NHOUR': write every NHOURSAVE
! 'END': write at end of run
! 'OUTDATE': write every OUTDATE(1:OUTDATE_NDUMP)
!
! To make a nested run:
! 1) run with MODE_SAVE='NHOUR' to write out 3d BC (name in filename_write defined below)
! 2) copy or link filename_write to filename_read_BC (for example "ln -s EMEP_OUT.nc EMEP_IN.nc")
! 3) run (in a smaller domain) with MODE_READ='NHOUR'
!
! Set MODE_SAVE/MODE_READ (in Nest_config nml) and out_DOMAIN (same namelist)
! Choose NHOURSAVE and NHOURREAD
! Also filename_read_BC and filename_read_3D should point to appropriate files
! Be careful to remove old BC files before making new ones.
!
! Grids may have any projection.
! Horizontal interpolation uses a weighted average of the four closest points
! This will work also if points in the present grid are not covered by the external grid.
! Vertical interpolation is done from hybrid coordinates.
!
!To do:
! It should be possible to save only xn_adv_bnd if the inner grid is known for the outer grid.
! The routines should be thought together with GlobalBC_mod (can it replace it?)
!----------------------------------------------------------------------------!
! External Boundary (BC) and Initial Conditions (IC)
! ExternalBICs_mod should handle different for different external sources.
! Experiment specific information must be set on ExternalBICs namelists.
! So far coded for CIFS (CAMS50/71) and EnsClimRCA(?) work.
use ExternalBICs_mod, only: set_extbic, icbc, ICBC_FMT,&
EXTERNAL_BIC_SET, EXTERNAL_BC, EXTERNAL_BIC_NAME, TOP_BC, &
iw, ie, js, jn, kt, & ! i West/East bnd; j North/South bnd; k Top
filename_eta,BC_DAYS
!----------------------------------------------------------------------------!
use CheckStop_mod, only: CheckStop,check=>CheckNC
use Chemfields_mod, only: xn_adv ! emep model concs.
use ChemDims_mod, only: NSPEC_ADV, NSPEC_SHL
use ChemSpecs_mod, only: species_adv
use GridValues_mod, only: A_mid,B_mid, glon,glat, i_fdom,j_fdom, RestrictDomain
use Io_mod, only: open_file,IO_TMP,IO_NML,PrintLog
use InterpolationRoutines_mod, only : grid2grid_coeff,point2grid_coeff
use MetFields_mod, only: roa
use Config_module, only: Pref,PT,KMAX_MID,MasterProc,NPROC,DataDir,GRID,&
IOU_INST,RUNDOMAIN,USES
use Debug_module, only: DEBUG_NEST,DEBUG_ICBC=>DEBUG_NEST_ICBC
use MPI_Groups_mod
use netcdf, only: nf90_open,nf90_write,nf90_close,nf90_inq_dimid,&
nf90_inquire_dimension,nf90_inq_varid,&
nf90_inquire_variable,nf90_get_var,nf90_get_att,&
nf90_put_att,nf90_noerr,nf90_nowrite,nf90_global
use netcdf_mod, only: Out_netCDF,&
CDFtype=>Real4,ReadTimeCDF,max_filename_length
use OwnDataTypes_mod, only: Deriv,TXTLEN_SHORT
use Par_mod, only: me,li0,li1,lj0,lj1,limax,ljmax,GIMAX,GJMAX,gi0,gj0,gi1,gj1
use Pollen_const_mod, only: pollen_check
use TimeDate_mod, only: date,current_date,nmdays
use TimeDate_ExtraUtil_mod, only: date2nctime,nctime2date,nctime2string,&
date2string,date2file,compare_date
use Units_mod, only: Units_Scale
use SmallUtils_mod, only: find_index,key2str,to_upper
use ChemGroups_mod, only: chemgroups
implicit none
! Nested input/output on OUTDATE mode
integer,private,parameter :: OUTDATE_NDUMP_MAX = 4 ! Number of nested output
integer, public, save :: OUTDATE_NDUMP = 1 ! Read by emepctm.f90
! on forecast run (1: start next forecast; 2-4: NMC statistics)
type(date), public :: outdate(OUTDATE_NDUMP_MAX)=date(-1,-1,-1,-1,-1)
!coordinates of subdomain to write, relative to FULL domain (only used in write mode)
integer, public, save :: out_DOMAIN(4) ! =[istart,iend,jstart,jend]
!/-- subroutines
public :: readxn
public :: wrtxn
private
logical, private, save :: mydebug = .false.
integer,private,parameter :: len_mode=20
character(len=len_mode),private, parameter :: &
READ_MODES(5)=[character(len=len_mode)::'NONE','RESTART','NHOUR','START','MONTH'],&
SAVE_MODES(5)=[character(len=len_mode)::'NONE','OUTDATE','NHOUR','END','MONTH']
character(len=len_mode),public, save :: &
MODE_READ='',& ! read mode
MODE_SAVE='' ! write mode
integer, private, save :: NHOURSAVE,NHOURREAD ! write/read frequency
!if(NHOURREAD<NHOURSAVE) the data is interpolated in time
character(len=max_filename_length),public, save :: &
template_read_3D = 'EMEP_IN.nc',& ! Different paths can be set here
template_read_BC = 'EMEP_IN.nc',& ! for each of the IO IC/BC files,
template_write = 'EMEP_OUT.nc' ! on Nest_config namelist, if needed.
character(len=max_filename_length),private, save :: &
filename_read_3D = 'template_read_3D',& ! Overwritten in readxn and wrtxn.
filename_read_BC = 'template_read_BC',& ! Filenames are updated according to date
filename_write = 'template_write' ! following respective templates
logical,private, save :: & ! if IC/BC are in the same model/run
native_grid_3D = .false.,& ! grid, the expensive call to
native_grid_BC = .false.,& ! grid2grid_coeff in init_nest can be avoided
omit_zero_write= .false. ! skip const=0.0 variables
character(len=max_filename_length),public, save :: MET_inner ='NOTSET' !path to metdata for inner grid
integer, save, public :: RUNDOMAIN_inner(4)=-1 ! RUNDOMAIN used for run in inner grid
! Limit output, e.g. for NMC statistics (3DVar)
character(len=TXTLEN_SHORT), private, save, dimension(NSPEC_ADV) :: &
WRITE_SPC = "" ! If these varables remain ""
character(len=TXTLEN_SHORT), private, save, dimension(size(chemgroups)) :: &
WRITE_GRP = "" ! all advected species will be written out.
real(kind=8), parameter :: &
halfsecond=0.5/(24.0*3600.0)! used to avoid rounding errors
!BC values at boundaries in present grid
real, save, allocatable, dimension(:,:,:,:) :: &
xn_adv_bndw, xn_adv_bnde, & ! west and east
xn_adv_bnds, xn_adv_bndn, & ! north and south
xn_adv_bndt ! top
real, save, allocatable, dimension(:) :: &
ndays_ext ! time stamps in days since 1900. NB: only defined on MasterProc
!dimension of external grid for BC
integer,save :: N_ext_BC ! Note: only defined on MasterProc
integer,save :: KMAX_ext_BC
integer,save :: itime
real(kind=8),save :: rtime_saved(2)
integer,save :: &
date_nextfile(4), & ! date corresponding to the next BC file to read
NHOURS_Stride_BC ! number of hours between start of two consecutive records in BC files
integer, public, parameter :: &
NHOURS_Stride_BC_default=6 !time between records if only one record per file (RCA for example)
type(icbc), private, target, dimension(NSPEC_ADV) :: &
adv_ic=icbc('none','none',1.0,.false.,.false.,-1) ! Initial 3D IC/BC: spcname,varname,wanted,found,ixadv
type(icbc), private, pointer, dimension(:) :: &
adv_bc=>null() ! Time dependent BC: spcname,varname,wanted,found,ixadv
logical, allocatable, save :: mask_restrict(:,:)
logical, save :: MASK_SET=.false.
contains
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
subroutine Config_Nest()
integer :: ios,i
logical, save :: first_call=.true.
NAMELIST /Nest_config/ MODE_READ,MODE_SAVE,NHOURREAD,NHOURSAVE, &
template_read_3D,template_read_BC,template_write,BC_DAYS,&
native_grid_3D,native_grid_BC,omit_zero_write,out_DOMAIN,&
MET_inner,RUNDOMAIN_inner,&
WRITE_SPC,WRITE_GRP,OUTDATE_NDUMP,outdate
if(.not.first_call)return
mydebug = DEBUG_NEST.and.MasterProc
! default write/read supported modes: do nothing
MODE_READ='NONE'
MODE_SAVE='NONE'
! write/read frequency: Hours between consecutive saves(wrtxn)/reads(readxn)
NHOURSAVE=3 ! Between wrtxn calls. Should be fraction of 24
NHOURREAD=1 ! Between readxn calls. Should be fraction of 24
! Default domain for write modes
out_DOMAIN=RUNDOMAIN
rewind(IO_NML)
read(IO_NML,NML=Nest_config,iostat=ios)
call CheckStop(ios,"NML=Nest_config")
if(mydebug)then
write(*,*) "NAMELIST IS "
write(*,NML=Nest_config)
end if
if(MODE_READ=='')then
MODE_READ='NONE'
else
MODE_READ=to_upper(MODE_READ)
end if
if(MODE_SAVE=='')then
MODE_SAVE='NONE'
else
MODE_SAVE=to_upper(MODE_SAVE)
end if
! write/read supported modes
if(MasterProc)then
call CheckStop(.not.any(MODE_READ==READ_MODES),&
"Config_Nest: Unsupported MODE_READ='"//trim(MODE_READ))
call CheckStop(.not.any(MODE_SAVE==SAVE_MODES),&
"Config_Nest: Unsupported MODE_SAVE='"//trim(MODE_SAVE))
end if
! write/read frequency should be fraction of 24
if(MasterProc)then
call CheckStop(mod(24,NHOURSAVE),"Config_Nest: NHOURSAVE should be fraction of 24")
call CheckStop(mod(24,NHOURREAD),"Config_Nest: NHOURREAD should be fraction of 24")
end if
! expand DataDir/GRID keyswords
template_read_3D=key2str(template_read_3D,'DataDir',DataDir)
template_read_3D=key2str(template_read_3D,'GRID',GRID)
template_read_BC=key2str(template_read_BC,'DataDir',DataDir)
template_read_BC=key2str(template_read_BC,'GRID',GRID)
template_write =key2str(template_write ,'DataDir',DataDir)
template_write =key2str(template_write ,'GRID',GRID)
! Update filenames according to date following templates defined on Nest_config
call init_icbc(cdate=current_date)
! Ensure out-domain is not larger than run-domain
if(MODE_SAVE/='NONE')call RestrictDomain(out_DOMAIN)
! Ensure that only OUTDATE_NDUMP are taking into account
if(MODE_SAVE=='OUTDATE')then
if(OUTDATE_NDUMP<OUTDATE_NDUMP_MAX)&
outdate(OUTDATE_NDUMP+1:OUTDATE_NDUMP_MAX)%day=0
if(MasterProc)&
write (*,"(1X,A,10(1X,A,:,','))")'OUTDATE nest/dump at:',&
(date2string("YYYY-MM-DD hh:mm:ss",outdate(i)),i=1,OUTDATE_NDUMP)
end if
first_call=.false.
end subroutine Config_Nest
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
subroutine readxn(indate)
type(date), intent(in) :: indate ! Gives year..seconds
integer :: n,i,j,k,KMAX_BC,bc,ndate(4)
real(kind=8):: ndays_indate
! real , dimension(48,48,20) ::data
real :: W1,W2
logical, save :: first_call=.true.
logical :: fexist_3D=.false.,fexist_BC=.false.
integer, save :: oldmonth=0
call Config_Nest()
if(mydebug) write(*,*)'Nest:Read BC, MODE=',MODE_READ
if(MODE_READ=='NONE')return
KMAX_BC=KMAX_MID
ndate(1:4)=[indate%year,indate%month,indate%day,indate%hour]
call date2nctime(ndate,ndays_indate)
if(first_call)date_nextfile=ndate
select case(MODE_READ)
case('MONTH') ! monthly input file
if(indate%month==oldmonth)return
if(MasterProc.and.oldmonth==0) write(*,*)'Nest: Initialzing IC'
oldmonth=indate%month
if(MasterProc) write(*,*)'Nest: New month, reset BC'
case('START')
if(.not.first_call)return
first_call=.false.
filename_read_3D=date2string(template_read_3D,ndate,mode='YMDH',debug=mydebug)
if(MasterProc) write(*,*)'Nest RESET ALL XN 3D ',trim(filename_read_3D)
call reset_3D(ndays_indate)
return
case default
!if(MasterProc) print *,'call to READXN',indate%hour,indate%seconds
if(mod(indate%hour,NHOURREAD)/=0.or.indate%seconds/=0)return
end select
! never comes to this point if MODE=100, 11 or 12
if(DEBUG_NEST.and.MasterProc) write(*,*) 'Nest: kt', kt, first_call
! Update filenames according to date following templates defined on Nest_config nml
filename_read_3D=date2string(template_read_3D,ndate,&
mode='YMDH',debug=mydebug)
filename_read_BC=date2file (template_read_BC,ndate,BC_DAYS,"days",&
mode='YMDH',debug=mydebug)
inquire(file=filename_read_3D,exist=fexist_3D)
inquire(file=filename_read_BC,exist=fexist_BC)
if(first_call)then
first_call=.false.
if(fexist_3D)then
if(MasterProc)write(*,*)'Nest RESET ALL XN 3D ',trim(filename_read_3D)
call reset_3D(ndays_indate)
else
if(MasterProc)write(*,*)'No Nest IC file found: ',trim(filename_read_3D)
end if
! the first hour only these values are used, no real interpolation between two records
if(fexist_BC)then
if(mydebug) write(*,*)'Nest: READING FIRST BC DATA from ',&
trim(filename_read_BC), ndays_indate
call read_newdata_LATERAL(ndays_indate)
if(mydebug) write(*,"(a,5i4)")'Nest: iw, ie, js, jn, kt ',iw,ie,js,jn,kt
end if
end if
if(.not.fexist_BC)then
if(MasterProc)write(*,*)'No Nest BC file found: ',trim(filename_read_BC)
return
end if
if(ndays_indate-rtime_saved(2)>halfsecond.or.MODE_READ=='MONTH')then
! look for a new data set
if(MasterProc) write(*,*)'Nest: READING NEW BC DATA from ',&
trim(filename_read_BC)
call read_newdata_LATERAL(ndays_indate)
end if
! make weights for time interpolation
if(MODE_READ=='MONTH')then ! don't interpolate for now
W1=0.0; W2=1.0 ! use last read value
else
W1=1.0; W2=0.0 ! default: use first read value
if(rtime_saved(2)-rtime_saved(1)<halfsecond)then
W1=0.0; W2=1.0 ! use last read value
elseif(ndays_indate-rtime_saved(1)>halfsecond)then
W2=(ndays_indate-rtime_saved(1))/(rtime_saved(2)-rtime_saved(1))
W1=1.0-W2 ! interpolate
end if
end if
if(DEBUG_NEST.and.MasterProc) then
write(*,*) 'Nesting BC 2D: time weights : ',W1,W2
write(*,*) 'Nesting BC 2D: time stamps : ',rtime_saved(1),rtime_saved(2)
end if
do bc=1,size(adv_bc)
if(.not.(adv_bc(bc)%wanted.and.adv_bc(bc)%found))cycle
n=adv_bc(bc)%ixadv
if(DEBUG_ICBC.and.MasterProc) write(*,"(2(A,1X),I0,'-->',I0)") &
'NestICBC: Nesting component',trim(adv_bc(bc)%varname),bc,n
forall (i=iw:iw, k=1:KMAX_BC, j=1:ljmax, i>=1) &
xn_adv(n,i,j,k)=W1*xn_adv_bndw(n,j,k,1)+W2*xn_adv_bndw(n,j,k,2)
forall (i=ie:ie, k=1:KMAX_BC, j=1:ljmax, i<=limax) &
xn_adv(n,i,j,k)=W1*xn_adv_bnde(n,j,k,1)+W2*xn_adv_bnde(n,j,k,2)
forall (j=js:js, k=1:KMAX_BC, i=1:limax, j>=1) &
xn_adv(n,i,j,k)=W1*xn_adv_bnds(n,i,k,1)+W2*xn_adv_bnds(n,i,k,2)
forall (j=jn:jn, k=1:KMAX_BC, i=1:limax, j<=ljmax) &
xn_adv(n,i,j,k)=W1*xn_adv_bndn(n,i,k,1)+W2*xn_adv_bndn(n,i,k,2)
forall (k=kt:kt, i=1:limax, j=1:ljmax, k>=1) &
xn_adv(n,i,j,k)=W1*xn_adv_bndt(n,i,j,1)+W2*xn_adv_bndt(n,i,j,2)
end do
call CheckStop(EXTERNAL_BIC_NAME=="RCA",&
"WORK NEEDED: RCA BICs commented out in Nest_mod - not consistent with all chem schemes")
end subroutine readxn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
subroutine wrtxn(indate,WriteNow)
type(date), intent(in) :: indate
logical, intent(in) :: WriteNow !Do not check indate value
real :: data(LIMAX,LJMAX,KMAX_MID) ! Data array
logical, parameter :: APPEND=.false.
type(Deriv) :: def1 ! definition of fields
integer :: n,i,j,k,iotyp,ndim,kmax,ncfileID
real :: scale
logical :: fexist, wanted, wanted_out, overwrite
logical, save :: first_call=.true.
call Config_Nest()
if(MODE_SAVE=='NONE')return
! Check if the file exist already at start of run. Do not wait until first write to stop!
! If you know what you are doing you can set paramter APPEND=.true.,
! and the new data will be appended to the file
overwrite=first_call.and..not.APPEND
if(overwrite.and.MasterProc)then
filename_write=date2string(template_write,indate,mode='YMDH',debug=mydebug)
inquire(file=fileName_write,exist=overwrite)
call CheckStop(overwrite.and.MODE_SAVE/='OUTDATE',&
"Nest: Refuse to overwrite. Remove this file: "//trim(fileName_write))
end if
select case(MODE_SAVE)
case('END')
if(.not.WriteNow)return
case('OUTDATE')
outdate(:)%seconds=0 ! output only at full hours
if(.not.compare_date(OUTDATE_NDUMP,indate,outdate(:OUTDATE_NDUMP),&
wildcard=-1))return
if(MasterProc) write(*,*)&
date2string(" Forecast nest/dump at YYYY-MM-DD hh:mm:ss",indate)
case('MONTH')
if(indate%month==1.or.indate%day/=1.or.indate%hour/=0.or.indate%seconds/=0)return
case default
if(mod(indate%hour,NHOURSAVE)/=0.or.indate%seconds/=0)return
end select
iotyp=IOU_INST
ndim=3 !3-dimensional
kmax=KMAX_MID
scale=1.0
def1%class='Advected' ! written
def1%avg=.false. ! not used
def1%index=0 ! not used
def1%scale=scale ! not used
def1%iotype='' ! not used
def1%name='' ! written
def1%unit='mix_ratio' ! written
! Update filenames according to date following templates defined on Nest_config nml
! e.g. set template_write="EMEP_BC_MMYYYY.nc" on namelist for different names each month
filename_write=date2string(template_write,indate,mode='YMDH',debug=mydebug)
if(MasterProc)then
inquire(file=fileName_write,exist=fexist)
write(*,*)'Nest:write data ',trim(fileName_write)
end if
CALL MPI_BCAST(fexist,1,MPI_LOGICAL,0,MPI_COMM_CALC,IERROR)
overwrite=fexist.and.first_call.and..not.APPEND
! Limit output, e.g. for NMC statistics (3DVar and restriction to inner grid BC)
if(first_call)then
call init_icbc(cdate=indate)
if(any([WRITE_GRP,WRITE_SPC]/=""))then
adv_ic(:)%wanted=.false.
do n=1,size(WRITE_GRP)
if(WRITE_GRP(n)=="")cycle
i=find_index(WRITE_GRP(n),chemgroups(:)%name)
if(i>0)then
where(chemgroups(i)%specs>NSPEC_SHL) &
adv_ic(chemgroups(i)%specs-NSPEC_SHL)%wanted=.true.
elseif(MasterProc)then
write(*,"(A,':',/2(2X,A,1X,'''',A,'''',1X,A,'.'))")&
"Warning (wrtxn)", &
"Wanted group",trim(WRITE_GRP(n)),"was not found", &
"Can not be written to file:",trim(filename_write),""
end if
end do
do n=1,size(WRITE_SPC)
if(WRITE_SPC(n)=="")cycle
i=find_index(WRITE_SPC(n),species_adv(:)%name)
if(i>0)then
adv_ic(i)%wanted=.true.
elseif((DEBUG_NEST.or.DEBUG_ICBC).and.MasterProc)then
write(*,"(A,':',/2(2X,A,1X,'''',A,'''',1X,A,'.'))")&
"Warning (wrtxn)", &
"Wanted specie",trim(WRITE_SPC(n)),"was not found", &
"Can not be written to file:",trim(filename_write),""
end if
end do
elseif(USES%POLLEN)then
! POLLEN group members are written to pollen restart/dump file
call pollen_check(igrp=i)
if(i>0)then
where(chemgroups(i)%specs>NSPEC_SHL) &
adv_ic(chemgroups(i)%specs-NSPEC_SHL)%wanted=.false.
if((DEBUG_NEST.or.DEBUG_ICBC).and.MasterProc)&
write(*,"(A,':',/2(2X,A,1X,'''',A,'''',1X,A,'.'))")&
"Warning (wrtxn)", &
"Group","POLLEN","is written to pollen restart/dump file", &
"Will not be written to file:",trim(filename_write),""
end if
end if
do n=1,NSPEC_ADV
if(.not.adv_ic(n)%wanted)then
if((DEBUG_NEST.or.DEBUG_ICBC).and.MasterProc)&
write(*,"(A,':',/2(2X,A,1X,'''',A,'''',1X,A,'.'))")&
"Nest(wrtxn) DEBUG_ICBC",&
"Variable",trim(species_adv(n)%name),"is not wanted as IC",&
"Will not be written to file:",trim(filename_write),""
elseif(omit_zero_write)then ! further reduce output
wanted=any(xn_adv(n,:,:,:)/=0.0)
CALL MPI_ALLREDUCE(wanted,wanted_out,1,MPI_LOGICAL,MPI_LOR,&
MPI_COMM_CALC,IERROR)
adv_ic(n)%wanted=wanted_out
if(.not.adv_ic(n)%wanted.and.&
(DEBUG_NEST.or.DEBUG_ICBC).and.MasterProc)&
write(*,"(A,':',/2(2X,A,1X,'''',A,'''',1X,A,'.'))")&
"Nest(wrtxn) DEBUG_ICBC",&
"Variable",trim(species_adv(n)%name),"was found constant=0.0",&
"Will not be written to file:",trim(filename_write),""
end if
end do
if(MET_inner /= "NOTSET")then
! find region that is really needed, i.e. boundaries of inner grid
!find lon and lat of inner grid restricted to BC
call init_mask_restrict(MET_inner,RUNDOMAIN_inner)
endif
end if
!do first one loop to define the fields, without writing them (for performance purposes)
ncfileID=-1 ! must be <0 as initial value
if(.not.fexist.or.overwrite)then
do n=1,NSPEC_ADV
if(.not.adv_ic(n)%wanted)cycle
def1%name=species_adv(n)%name ! written
!! data=xn_adv(n,:,:,:)
call Out_netCDF(iotyp,def1,ndim,kmax,data,scale,CDFtype=CDFtype,&
out_DOMAIN=out_DOMAIN,create_var_only=.true.,overwrite=overwrite,&
fileName_given=trim(fileName_write),ncFileID_given=ncFileID)
overwrite=.false.
end do
end if
do n=1,NSPEC_ADV
if(.not.adv_ic(n)%wanted)cycle
def1%name=species_adv(n)%name ! written
if(MASK_SET)then
do k=1,KMAX_MID
do j=1,LJMAX
do i=1,LIMAX
if(mask_restrict(i,j))then
data(i,j,k)=xn_adv(n,i,j,k)
else
data(i,j,k)=0.0
endif
end do
end do
end do
else
data=xn_adv(n,:,:,:)
endif
call Out_netCDF(iotyp,def1,ndim,kmax,data,scale,CDFtype=CDFtype,&
out_DOMAIN=out_DOMAIN,create_var_only=.false.,&
fileName_given=trim(fileName_write),ncFileID_given=ncFileID)
end do
if(first_call .and. MET_inner /= "NOTSET" .and. me==0)then
!mark the file as defined in a restricted area only
call check(nf90_put_att(ncFileID,nf90_global,"restricted","BC_restricted"))
endif
first_call=.false.
if(MasterProc)call check(nf90_close(ncFileID))
end subroutine wrtxn
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
subroutine init_icbc(idate,cdate,ndays,nsecs)
!----------------------------------------------------------------------------!
! Setup IC/BC detailed description.
! ICs are assumed to come from emepctm.
!
! adv_ic IC detailed description for all adv species
! adv_bc BC detailed description relevant adv species
! EXTERNAL_BC External (non emepctm) BC detailed description/setup
! EXTERNAL_BIC_SET EXTERNAL_BC has been set (adv_bc=>EXTERNAL_BC)
! otherwise Assume emepctm BCs (adv_bc=>adv_ic)
!----------------------------------------------------------------------------!
integer, intent(in), optional :: idate(4)
type(date),intent(in), optional :: cdate
real(kind=8),intent(in),optional:: ndays
integer, intent(in),optional:: nsecs
logical, save :: first_call=.true.
integer :: n,dat(4)
if(.not.first_call)return
first_call=.false.
! One of the date formats needs to be provided
call CheckStop(count([present(idate),present(cdate),present(ndays),&
present(nsecs)]),1,"init_icbc: wrong date option")
! Update filenames according to date following templates defined on Nest_config nml
if(present(idate)) dat=idate
if(present(cdate)) dat=[cdate%year,cdate%month,cdate%day,cdate%hour]
if(present(ndays)) call nctime2date(dat,ndays)
if(present(nsecs)) call nctime2date(dat,nsecs)
call set_extbic(dat) ! set mapping, EXTERNAL_BC, TOP_BC
if(.not.EXTERNAL_BIC_SET.and.MODE_READ=='NONE'.and.MODE_SAVE=='NONE')return !No nesting
filename_read_3D=date2string(template_read_3D,dat,&
mode='YMDH',debug=mydebug)
filename_read_BC=date2file (template_read_BC,dat,BC_DAYS,"days",&
mode='YMDH',debug=mydebug)
filename_write =date2string(template_write ,dat,&
mode='YMDH',debug=mydebug)
adv_ic(:)%ixadv=(/(n,n=1,NSPEC_ADV)/)
adv_ic(:)%spcname=species_adv(:)%name
adv_ic(:)%varname=species_adv(:)%name
adv_ic(:)%frac=1.0
adv_ic(:)%wanted=.true.
adv_ic(:)%found=find_icbc(filename_read_3D,adv_ic%varname(:))
if(EXTERNAL_BIC_SET) then
adv_bc=>EXTERNAL_BC
adv_bc(:)%found=find_icbc(filename_read_bc,adv_bc%varname(:))
else
adv_bc=>adv_ic
end if
if(MasterProc)then
do n = 1,size(adv_ic%varname)
if(.not.adv_ic(n)%found)then
call PrintLog("WARNING: IC variable '"//trim(adv_ic(n)%varname)//"' not found")
elseif(DEBUG_NEST.or.DEBUG_ICBC)then
write(*,*) "init_icbc filled adv_ic "//trim(adv_ic(n)%varname)
end if
end do
do n = 1,size(adv_bc%varname)
if(.not.adv_bc(n)%found)then
call PrintLog("WARNING: BC variable '"//trim(adv_bc(n)%varname)//"' not found")
elseif(DEBUG_NEST.or.DEBUG_ICBC)then
write(*,*) "init_icbc filled adv_bc "//trim(adv_bc(n)%varname)
end if
end do
end if
if((DEBUG_NEST.or.DEBUG_ICBC).and.MasterProc)then
write(*,"(a)") "Nest: DEBUG_ICBC Variables:",&
trim(filename_read_3D),trim(filename_read_BC)
write(*,"((1X,A,I3,'->',"//ICBC_FMT//"))") &
('Nest: ADV_IC',n,adv_ic(n),n=1,size(adv_ic)),&
('Nest: ADV_BC',n,adv_bc(n),n=1,size(adv_bc))
end if
contains
function find_icbc(filename_read,varname) result(found)
!----------------------------------------------------------------------------!
! Check if variables (varname) are present on file (filename_read)
!----------------------------------------------------------------------------!
implicit none
character(len=*), intent(in) :: filename_read
character(len=*), dimension(:), intent(in) :: varname
logical, dimension(size(varname)) :: found
integer :: ncFileID,varID,status,n
found(:)=.false.
if(MasterProc)then
status=nf90_open(trim(filename_read),nf90_nowrite,ncFileID)
if(status/=nf90_noerr) then
print *,'icbc: not found ',trim(filename_read)
else
print *,'icbc: reading ',trim(filename_read)
do n=1,size(varname)
if(varname(n)=="") cycle
status=nf90_inq_varid(ncFileID,trim(varname(n)),varID)
found(n)=(status==nf90_noerr)
end do
call check(nf90_close(ncFileID))
end if
end if
CALL MPI_BCAST(found,size(found),MPI_LOGICAL,0,MPI_COMM_CALC,IERROR)
end function find_icbc
end subroutine init_icbc
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
subroutine init_nest(ndays_indate,filename_read,native_grid,IIij,JJij,Weight,&
k1_ext,k2_ext,weight_k1,weight_k2,N_ext,KMAX_ext,GIMAX_ext,GJMAX_ext)
logical, parameter :: USE_LAST_HYBRID_LEVELS=.true.
character(len=*),intent(in) :: filename_read
logical,intent(in) :: native_grid
real ,intent(out):: Weight(4,LIMAX,LJMAX)
integer ,intent(out)::IIij(4,LIMAX,LJMAX),JJij(4,LIMAX,LJMAX)
integer, intent(out), dimension(*) :: k1_ext,k2_ext
real, intent(out), dimension(*) :: weight_k1,weight_k2
integer ,intent(out)::N_ext,KMAX_ext,GIMAX_ext,GJMAX_ext
real(kind=8) :: ndays_indate
integer :: ncFileID,timeDimID,varid,status,dimIDs(3)
real :: P_emep
integer :: i,j,k,n,k_ext
real, allocatable, dimension(:,:) ::lon_ext,lat_ext
real, allocatable, dimension(:) ::hyam,hybm,P_ext,temp_ll
character(len=80) ::projection,word,iDName,jDName
logical :: reversed_k_BC,time_exists,fexist
rtime_saved = -99999.9 !initialization
!Read dimensions (global)
if(MasterProc)then
status = nf90_open(trim(filename_read),nf90_nowrite,ncFileID)
if(status/=nf90_noerr) then
print *,'init_Nest: not found',trim(filename_read)
return
else
print *,'init_Nest: reading ',trim(filename_read)
end if
projection='Unknown'
status = nf90_get_att(ncFileID,nf90_global,"projection",projection)
if(status==nf90_noerr) then
if(projection=='lon_lat')projection='lon lat'
write(*,*)'Nest: projection: '//trim(projection)
else
projection='lon lat'
write(*,*)'Nest: projection not found for ',&
trim(filename_read)//', assuming '//trim(projection)
end if
!get dimensions id/name/len: include more dimension names, if necessary
GIMAX_ext=get_dimLen([character(len=12)::"i","lon","longitude"],name=iDName)
GJMAX_ext=get_dimLen([character(len=12)::"j","lat","latitude" ],name=jDName)
KMAX_ext =get_dimLen([character(len=12)::"k","mlev","lev","level"])
select case(projection)
case('Stereographic')
call CheckStop("i",iDName,"Nest: unsuported "//&
trim(iDName)//" as i-dimension on "//trim(projection)//" projection")
call CheckStop("j",jDName,"Nest: unsuported "//&
trim(jDName)//" as j-dimension on "//trim(projection)//" projection")
case('lon lat')
call CheckStop("lon",iDName(1:3),"Nest: unsuported "//&
trim(iDName)//" as i-dimension on "//trim(projection)//" projection")
call CheckStop("lat",jDName(1:3),"Nest: unsuported "//&
trim(jDName)//" as j-dimension on "//trim(projection)//" projection")
case default
!call CheckStop("Nest: unsuported projection "//trim(projection))
!write(*,*)'GENERAL PROJECTION ',trim(projection)
call CheckStop("i",iDName,"Nest: unsuported "//&
trim(iDName)//" as i-dimension on "//trim(projection)//" projection")
call CheckStop("j",jDName,"Nest: unsuported "//&
trim(jDName)//" as j-dimension on "//trim(projection)//" projection")
end select
N_ext=0
status = nf90_inq_dimid(ncFileID,"time",timeDimID)
time_exists=(status==nf90_noerr)
if(time_exists) then
call check(nf90_inquire_dimension(ncFileID,timedimID,len=N_ext))
else
status = nf90_inq_dimid(ncFileID,"Months",dimID=timeDimID)
if(status==nf90_noerr)then
call check(nf90_inquire_dimension(ncFileID,timedimID,len=N_ext))
call CheckStop(N_ext,12,'Nest BC: did not find 12 months')
else
write(*,*)'Nest: time dimension not found. Assuming only one record '
N_ext=1
end if
end if
write(*,*)'Nest: dimensions external grid',GIMAX_ext,GJMAX_ext,KMAX_ext,N_ext
if(.not.allocated(ndays_ext))then
allocate(ndays_ext(N_ext))
elseif(size(ndays_ext)<N_ext)then
if(Masterproc)write(*,*)'Nest: Sizes times ',N_ext
deallocate(ndays_ext)
allocate(ndays_ext(N_ext))
end if
end if
CALL MPI_BCAST(GIMAX_ext,4*1,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
CALL MPI_BCAST(GJMAX_ext,4*1,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
CALL MPI_BCAST(KMAX_ext ,4*1,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
!CALL MPI_BCAST(N_ext,4*1,MPI_BYTE,0,MPI_COMM_CALC,IERROR) !not needed by others than MasterProc
allocate(lon_ext(GIMAX_ext,GJMAX_ext),lat_ext(GIMAX_ext,GJMAX_ext))
allocate(hyam(KMAX_ext+1),hybm(KMAX_ext+1),P_ext(KMAX_ext))
if(MasterProc)then
!Read lon lat of the external grid (global)
if(trim(projection)==trim('lon lat')) then
call check(nf90_inq_varid(ncFileID,iDName,varID),&
"Read lon-variable: "//trim(iDName))
allocate(temp_ll(GIMAX_ext))
call check(nf90_get_var(ncFileID,varID,temp_ll))
lon_ext=SPREAD(temp_ll,2,GJMAX_ext)
deallocate(temp_ll)
call check(nf90_inq_varid(ncFileID,jDName,varID),&
"Read lat-variable: "//trim(jDName))
allocate(temp_ll(GJMAX_ext))
call check(nf90_get_var(ncFileID,varID,temp_ll))
lat_ext=SPREAD(temp_ll,1,GIMAX_ext)
deallocate(temp_ll)
else
call check(nf90_inq_varid(ncFileID,"lon",varID),"dim:lon")
call check(nf90_get_var(ncFileID,varID,lon_ext),"get:lon")
call check(nf90_inq_varid(ncFileID,"lat",varID),"dim:lat")
call check(nf90_get_var(ncFileID,varID,lat_ext),"get:lat")
end if
if(time_exists)then
call ReadTimeCDF(filename_read,ndays_ext,N_ext)
else
!cannot read time on file. assumes it is correct
ndays_ext(1)=ndays_indate
end if
if(MODE_READ=='MONTH')then
!asuming 12 monthes for BC, and 12 or 1 values for IC
ndays_ext(1)=0
do n=2,N_ext
ndays_ext(n)=ndays_ext(n-1)+nmdays(n-1)
end do
elseif(ndays_ext(1)-ndays_indate>halfsecond)then
write(*,*)'WARNING: Nest did not find BIC for date ',&
nctime2string('YYYY-MM-DD hh:mm:ss',ndays_indate)
write(*,*)'Nest first date found ',&
nctime2string('YYYY-MM-DD hh:mm:ss',ndays_ext(1))
end if
if(N_ext>1)then
NHOURS_Stride_BC = nint((ndays_ext(2)-ndays_ext(1))*24)
else
!use manually set stride:
NHOURS_Stride_BC = NHOURS_Stride_BC_default
end if
write(*,*)'Nest: new BC record every ',NHOURS_Stride_BC,' hours'
! Read pressure for vertical levels
write(*,*)'Nest: reading vertical levels'
status = nf90_inq_varid(ncFileID,"hyam",varID)
if(status==nf90_noerr) then
write(*,*)'Found hyam type levels (values at level midpoints)'
call check(nf90_inquire_variable(ncFileID,varID,dimIDs=dimIDs))
call check(nf90_inquire_dimension(ncFileID,dimIDs(1),len=k))
call CheckStop(k<KMAX_ext,"Nest BC, wrong hyam/hybm dimension")
if(USE_LAST_HYBRID_LEVELS)then
k_ext=1+k-KMAX_ext ! for 1+k-KMAX_ext .. k levels
else
k_ext=1 ! for 1 .. KMAX_ext levels
end if
if(k/=KMAX_ext)&
write(*,"(A,4(1X,A,I0))")'Nest BC warning:',&
'kdim #lev=',KMAX_ext,'and hyam/hybm #lev=',k,&
'. Using only levels ',k_ext,'..',k
call check(nf90_get_var(ncFileID,varID,hyam,start=(/k_ext/),count=(/KMAX_ext/)))
status = nf90_get_att(ncFileID,VarID,"units",word)
if(status==nf90_noerr)then
if(word(1:3)=='hPa')then
write(*,*)'Changing hyam from hPa to Pa'
hyam=100*hyam
end if
end if
call check(nf90_inq_varid(ncFileID,"hybm",varID))
call check(nf90_get_var(ncFileID,varID,hybm,start=(/k_ext/),count=(/KMAX_ext/)))
else
status = nf90_inq_varid(ncFileID,"hyai",varID)
if(status==nf90_noerr) then
write(*,*)'Found hyai type levels (values at level interfaces)'
call check(nf90_inquire_variable(ncFileID,varID,dimIDs=dimIDs))
call check(nf90_inquire_dimension(ncFileID,dimIDs(1),len=k))
call CheckStop(k<KMAX_ext+1,"Nest BC, wrong hyai/hybi dimension")
if(USE_LAST_HYBRID_LEVELS)then
k_ext=1+(k-1)-KMAX_ext ! for 1+k-KMAX_ext .. k levels
else
k_ext=1 ! for 1 .. KMAX_ext levels
end if
if(k/=KMAX_ext+1.and.MasterProc)&
write(*,"(A,4(1X,A,I0))")'Nest BC warning:',&
'kdim #lev=',KMAX_ext,'and hyam/hybm #lev=',k,&
'. Using only levels ',k_ext,'..',k
call check(nf90_get_var(ncFileID,varID,hyam,start=(/k_ext/),count=(/KMAX_ext+1/)))
status = nf90_get_att(ncFileID,VarID,"units",word)
if(status==nf90_noerr)then
if(word(1:3)=='hPa')then
write(*,*)'Changing hyai from hPa to Pa'
hyam=100*hyam
end if
end if
call check(nf90_inq_varid(ncFileID,"hybi",varID))
call check(nf90_get_var(ncFileID,varID,hybm,start=(/k_ext/),count=(/KMAX_ext+1/)))
do k=1,KMAX_ext
hyam(k)=0.5*(hyam(k)+hyam(k+1))
hybm(k)=0.5*(hybm(k)+hybm(k+1))
end do
else
inquire(file=filename_eta,exist=fexist)
status = nf90_inq_varid(ncFileID,"k",varID)
if(status==nf90_noerr) then
write(*,*)'Nest: assuming sigma level and PT=',PT,KMAX_ext
call check(nf90_get_var(ncFileID, varID, hybm,count=(/ KMAX_ext /) ))!NB: here assume = sigma
do k=1,KMAX_ext
hyam(k)=PT*(1.0-hybm(k))
end do
elseif(fexist) then
!read eta levels from ad-hoc text file
write(*,*)'Nest: Reading vertical level from ',trim(filename_eta)
call open_file(IO_TMP,"r",trim(filename_eta),needed=.true.)
do n=1,10000
read(IO_TMP,*)word
if(trim(word)=='vct')exit
end do
read(IO_TMP,*)(hyam(k),k=1,KMAX_ext+1)!NB: here = A_bnd, not mid
read(IO_TMP,*)(hybm(k),k=1,KMAX_ext+1)!NB: here = B_bnd, not mid
close(IO_TMP)
!convert to mid levels coefficients
do k=1,KMAX_ext
hyam(k)=0.5*(hyam(k)+hyam(k+1))
hybm(k)=0.5*(hybm(k)+hybm(k+1))
end do
else
status = nf90_inq_varid(ncFileID,"lev",varID)
if(status == nf90_noerr) then
call CheckStop('Pressure levels not yet implemented')
write(*,*)'Nest: assuming pressure levels and hPa'
call check(nf90_get_var(ncFileID, varID, hyam,count=(/ KMAX_ext /) ))
hyam=100.0*hyam ! hPa ->Pa
hybm=0.0
else
call CheckStop('Vertical coordinate Unknown/Not yet implemented')
end if
end if
end if
end if
call check(nf90_close(ncFileID))
end if !end MasterProc
CALL MPI_BCAST(lon_ext,8*GIMAX_ext*GJMAX_ext,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
CALL MPI_BCAST(lat_ext,8*GIMAX_ext*GJMAX_ext,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
CALL MPI_BCAST(hyam ,8*KMAX_ext ,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
CALL MPI_BCAST(hybm ,8*KMAX_ext ,MPI_BYTE,0,MPI_COMM_CALC,IERROR)
! find horizontal interpolation constants
! note that i,j are local and but IIij,JJij refer to the full nest-file
if(native_grid)then ! nest-file is on the model/run grid
forall(i=1:limax,j=1:ljmax)
IIij(:,i,j)=i_fdom(i)-RUNDOMAIN(1)+1
JJij(:,i,j)=j_fdom(j)-RUNDOMAIN(3)+1
Weight(:,i,j)=[1.0,0.0,0.0,0.0]
endforall
i=IIij(1,limax,ljmax);j=JJij(1,limax,ljmax)
call CheckStop((i>GIMAX_ext).or.(j>GJMAX_ext),&
'Nest: domain mismatch for native_grid')
else ! find the four closest points
call grid2grid_coeff(glon,glat,IIij,JJij,Weight,lon_ext,lat_ext,&
GIMAX_ext,GJMAX_ext,LIMAX,LJMAX,limax,ljmax,mydebug,1,1)
end if
deallocate(lon_ext,lat_ext)
!find vertical interpolation coefficients
!use pressure as reference
!we want, if possible, P_ext(k1) and P_ext(k2) to be on each side of P_emep
!We assume constant surface pressure, both for emep and external grid; should not be so
! important as long as they are both terrain following.
do k_ext=1,KMAX_EXT
P_ext(k_ext)=hyam(k_ext)+hybm(k_ext)*Pref
if(mydebug) write(*,fmt="(A,I3,F10.2)")'Nest: P_ext',k_ext,P_ext(k_ext)
end do
reversed_k_BC=(P_ext(1)>P_ext(2))
! .true. --> assumes k_ext=KMAX_EXT is top and k_ext=1 is surface
! .false. --> assumes k_ext=1 is top and k_ext=KMAX_EXT is surface
if(reversed_k_BC)then
do k=1,KMAX_MID
P_emep=A_mid(k)+B_mid(k)*Pref !Pa
if(mydebug) write(*,fmt="(A,I3,F10.2)")'Nest: P_emep',k,P_emep
!largest available P smaller than P_emep (if possible)
k1_ext(k)=1 !start at surface, and go up until P_emep
do k_ext=1,KMAX_EXT
if(P_ext(k_ext)<P_emep)exit
k1_ext(k)=k_ext
end do
!smallest available P larger than P_emep (if possible)
k2_ext(k)=KMAX_EXT !start at top, and go down until P_emep
if(k2_ext(k)==k1_ext(k))k2_ext(k)=KMAX_EXT-1 !avoid k2=k1
do k_ext=KMAX_EXT,1,-1
if(P_ext(k_ext)>P_emep)exit
if(k_ext/=k1_ext(k))k2_ext(k)=k_ext
end do
weight_k1(k)=(P_emep-P_ext(k2_ext(k)))/(P_ext(k1_ext(k))-P_ext(k2_ext(k)))
weight_k2(k)=1.0-weight_k1(k)
if(mydebug)&
write(*,fmt="(A,I3,2(A,I2,A,F5.2))")'Nest: level',k,&
' is the sum of level ',k1_ext(k),' weight ',weight_k1(k),&
' and level ',k2_ext(k),' weight ',weight_k2(k)
end do
else
do k=1,KMAX_MID
P_emep=A_mid(k)+B_mid(k)*Pref !Pa
if(mydebug) write(*,fmt="(A,I3,F10.2)")'Nest: P_emep',k,P_emep
!largest available P smaller than P_emep (if possible)
k1_ext(k)=KMAX_EXT !start at surface, and go up until P_emep
do k_ext=KMAX_EXT,1,-1
if(P_ext(k_ext)<P_emep)exit
k1_ext(k)=k_ext
end do
!smallest available P larger than P_emep (if possible)
k2_ext(k)=1 !start at top, and go down until P_emep
if(k2_ext(k)==k1_ext(k))k2_ext(k)=2 !avoid k2=k1
do k_ext=1,KMAX_EXT
if(P_ext(k_ext)>P_emep)exit
if(k_ext/=k1_ext(k))k2_ext(k)=k_ext
end do
weight_k1(k)=(P_emep-P_ext(k2_ext(k)))/(P_ext(k1_ext(k))-P_ext(k2_ext(k)))
weight_k2(k)=1.0-weight_k1(k)
if(mydebug) &
write(*,fmt="(A,I3,2(A,I2,A,F5.2))")'Nest: level',k,&
' is the sum of level ', k1_ext(k),' weight ',weight_k1(k),&
' and level ', k2_ext(k),' weight ',weight_k2(k)
end do
end if
deallocate(P_ext,hyam,hybm)
if(mydebug) &
write(*,*)'Nest: finished determination of interpolation parameters'
contains
function get_dimLen(dimName,id,name) result(len)
character(len=*), dimension(:), intent(in) :: dimName
integer, optional, intent(out):: id
character(len=*), optional, intent(out):: name
integer :: d, dID, len
do d=1,size(dimName)
status = nf90_inq_dimid(ncFileID,dimName(d),dID)
if(status==nf90_noerr)then
if(present(id)) id=dID
if(present(name))name=trim(dimName(d))
call check(nf90_inquire_dimension(ncFileID,dID,len=len),&
"get_dimLen: "//trim(dimName(d)))
exit
end if
end do
call CheckStop(status,nf90_noerr,'Nest: '//&
trim(dimName(1))//'-dimension not found: '//&
trim(filename_read)//'. Include new name in init_nest')