-
Notifications
You must be signed in to change notification settings - Fork 0
/
GridValues_mod.f90
executable file
·2900 lines (2585 loc) · 106 KB
/
GridValues_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 GridValues_mod
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
! Define parameters, variables and transformnations associated with grid and
! projection.
!
! Nomenclature:
! fulldomain is the largest grid, usually where metdata is defined.
! rundomain is a grid where the run is performed, smaller than fulldomain.
! subdomain: the domain covered by one MPI process or processor.
! restricted domain is a grid smaller than rundomain, where data is outputed;
! (the restricted domains are for instance, fullrun_DOMAIN,month_DOMAIN,
! day_DOMAIN,hour_DOMAIN).
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
use CheckStop_mod, only: CheckStop,StopAll,check=>CheckNC
use Functions_mod, only: great_circle_distance
use Io_Nums_mod, only: IO_LOG,IO_TMP
use MetFields_mod
use Config_module, only: &
KMAX_BND, KMAX_MID, & ! vertical extent
MasterProc,NPROC,IIFULLDOM,JJFULLDOM,RUNDOMAIN, JUMPOVER29FEB,&
PT,Pref,NMET,USE_EtaCOORDINATES,MANUAL_GRID,USE_WRF_MET_NAMES,&
startdate,NPROCX,NPROCY,Vertical_levelsFile,&
EUROPEAN_settings, GLOBAL_settings,USES,FORCE_PFT_MAPS_FALSE
use Debug_module, only: DEBUG ! -> DEBUG%GRIDVALUES
use MPI_Groups_mod!, only : MPI_BYTE, MPI_DOUBLE_PRECISION, MPI_LOGICAL, &
! MPI_MIN, MPI_MAX, &
! MPI_COMM_CALC, MPI_COMM_WORLD, MPISTATUS, IERROR, &
! ME_MPI, NPROC_MPI, ME_CALC, largeLIMAX,largeLJMAX
use Par_mod, only : &
GIMAX,GJMAX, & ! Size of rundomain
IRUNBEG,JRUNBEG, & ! start of rundomain in fulldomain coordinates
gi0,gj0, & ! rundomain coordinates of subdomain lower l.h. corner
gi1,gj1, & ! rundomain coordinates of subdomain uppet r.h. corner
limax,ljmax, & ! max i,j in this subdomain (can differ on other subdomains)
li0,li1,lj0,lj1, & ! start and end of i,j excluding outer frame of rundomain.
! li0=1 or 2, li1=limax or limax-1
me, & ! local processor
neighbor,WEST,EAST,SOUTH,NORTH,NOPROC, &
parinit,parinit_groups, &
MAXLIMAX,MAXLJMAX,MINLIMAX,MINLJMAX,tljmax,tlimax
use PhysicalConstants_mod, only: GRAV,PI,EARTH_RADIUS,deg2rad,rad2deg
use TimeDate_mod, only: current_date,date,Init_nmdays,nmdays
use TimeDate_ExtraUtil_mod, only: date2string
use InterpolationRoutines_mod, only: inside_1234
use netcdf, only: &
NF90_OPEN,NF90_NOWRITE,NF90_NOERR,NF90_CLOSE,&
NF90_GET_ATT,NF90_GLOBAL,NF90_INQ_DIMID,NF90_INQUIRE_DIMENSION,&
NF90_INQ_VARID,NF90_GET_VAR
implicit none
private
!-- contains subroutine:
public :: DefDebugProc ! => sets debug_proc, debug_li, debug_lj
public :: ij2lbm ! polar stereo grid to longitude latitude
public :: lb2ijm ! longitude latitude to grid in polar stereo
public :: ij2ijm ! polar grid1 to polar grid2
public :: lb2ij ! longitude latitude to (i,j) in any grid projection
public :: ij2lb ! polar stereo grid to longitude latitude
public :: lb_rot2lb !rotated lon lat to lon lat
public :: lb2UTM ! lon lat to Transverse Mercator
public :: UTM2lb ! Transverse Mercator to lon lat
interface lb2ij
module procedure lb2ij_real,lb2ij_int
end interface
private :: lb2ij_real,lb2ij_int
public :: coord_check ! normalize longitudes
public :: &
coord_in_gridbox, & ! Are coord (lon/lat) inside gridbox(i,j)?
coord_in_processor,& ! Are coord (lon/lat) inside local domain?
coord_in_domain ! Are coord (lon/lat) inside "domain" (full, run or sub)?
public :: RestrictDomain ! mask from full domain to rundomain
public :: GridRead
public :: extendarea_N ! returns array which includes neighbours from other subdomains
public :: set_EuropeanAndGlobal_Config
private :: Alloc_GridFields
private :: GetFullDomainSize
private :: find_poles
!** 1) Public (saved) Variables from module:
! Polar stereographic projection parameters
real, public, save :: &
xp=0.0, yp=1.0, & ! Coordinates of North pole
fi=0.0, & ! projections rotation angle around y axis
AN=1.0, & ! Distance on the map from pole to equator (No. of cells)
GRIDWIDTH_M=1.0, & ! width of grid at ref_latitude, in meters
ref_latitude =60. ! latitude at which projection is true (degrees)
! Rotated_Spherical grid prarameters
real, public, save :: &
grid_north_pole_latitude,grid_north_pole_longitude,&
dx_rot,dx_roti,x1_rot,y1_rot
! Lambert conformal projection parameters
real, public, save :: &
lon0_lambert,& ! reference longitude, also called phi, at which y=0 if lat=lat0_lambert
lat0_lambert,& ! reference latitude, at which x=0
lat_stand1_lambert,&! standard latitude at which mapping factor=1
lat_stand2_lambert,&! second standard latitude
y0_lambert,& ! reference y coordinate, also called rho0
k_lambert,& ! also called n, = sin(lat_stand1_lambert)
earth_radius_lambert,&! earth_radius used to define x and y in the met file. NOT USED
F_lambert,& ! normalization constant = cos(dr*lat0_lambert)*tan(PI/4+dr2*lat0_lambert)**k_lambert/k_lambert
x1_lambert,& ! x value at i=1
y1_lambert ! y value at j=1
!/ Variables to define full-domain (fdom) coordinates of local i,j values,
! and reciprocal variables.
integer, public, allocatable, save, dimension(:) :: &
i_fdom,j_fdom, & ! fdom coordinates of local i,j
i_local,j_local ! local coordinates of full-domain i,j
!Parameters for Vertical Hybrid coordinates:
real, public, save,allocatable, dimension(:) :: &
A_bnd,B_bnd,& ! first [Pa],second [1] constants at layer boundary
! (i.e. half levels in EC nomenclature)
A_bnd_met,B_bnd_met,& ! first [Pa],second [1] constants at layer boundary
! (i.e. half levels in EC nomenclature)
A_mid,B_mid,& ! first [Pa],second [1] constants at middle of layer
! (i.e. full levels in EC nomenclature)
dA,dB,& ! A_bnd(k+1)-A_bnd(k) [Pa],B_bnd(k+1)-B_bnd(k) [1]
! P = A + B*PS; eta = A/Pref + B
dEta_i,& ! 1/deta = 1/(dA/Pref + dB)
Eta_bnd,Eta_mid,& ! boundary,midpoint of eta layer
sigma_bnd,sigma_mid ! boundary,midpoint of sigma layer
real, public, save,allocatable, dimension(:,:) :: &
glon ,glat ,& !longitude,latitude of gridcell centers
gl_stagg ,gb_stagg,& !longitude,latitude of gridcell corners
!NB: gl_stagg,gb_stagg are here defined as the average of the four
! surrounding glat,glon.
! These differ slightly from the staggered points in the (i,j) grid.
rot_angle
real, public, save :: gbacmax,gbacmin,glacmax,glacmin
! EMEP grid definitions (old and official)
real, public, parameter :: &
xp_EMEP_official=8.,yp_EMEP_official=110.0,fi_EMEP=-32.0,&
ref_latitude_EMEP=60.0,GRIDWIDTH_M_EMEP=50000.0,&
an_EMEP=237.7316364, &! = 6.370e6*(1.0+0.5*sqrt(3.0))/50000.
xp_EMEP_old=43.0,yp_EMEP_old=121.0
!*** Map factor stuff:
real, public, save,allocatable, dimension(:,:) :: &
xm_i, & ! map-factor in i direction, between cell j and j+1
xm_j, & ! map-factor in j direction, between cell i and i+1
xm2, & ! xm*xm: area factor in the middle of a cell (i,j)
xmd, & ! 1/xm2
xm2ji,xmdji
!vertical "map factors"
real, public, save, allocatable,dimension(:) :: dhs1, dhs1i, dhs2i
!*** Grid Area
real, public, save,allocatable, dimension(:,:) :: GridArea_m2
integer, public, save :: &
debug_li=-99, debug_lj=-99 ! Local Coordinates of debug-site
logical, public, save :: debug_proc ! Processor with debug-site
character(len=100),public :: projection
integer, public, parameter :: MIN_ADVGRIDS = 5 !minimum size of a subdomain
integer, public :: Poles(2) ! Poles(1)=1 if North pole is found, Poles(2)=1:SP
integer, public :: Pole_Singular ! Pole_included=1 or 2 if the grid include
! at least one pole and has lat lon projection
logical, public :: Grid_Def_exist
integer, allocatable, save, public :: k1_met(:),k2_met(:)
real, allocatable, save, public :: x_k1_met(:)
logical, public, save :: External_Levels_Def=.false.
integer, public, save :: KMAX_MET !number of vertical levels from the meteo files
real, private :: u(2),v(2)!array for temporary use
contains
subroutine GridRead(meteo,cyclicgrid)
! the subroutine reads the grid parameters (projection, resolution etc.)
! defined by the meteorological fields
implicit none
character(len=*),intent(in):: meteo ! template for meteofile
integer, intent(out) :: cyclicgrid
integer :: nyear,nmonth,nday,nhour,k,ios
integer :: MIN_GRIDS
character(len=len(meteo)) :: filename !name of the input file
logical :: Use_Grid_Def=.false.!Experimental for now
nyear=startdate(1)
nmonth=startdate(2)
nday=startdate(3)
nhour=startdate(4)
current_date = date(nyear, nmonth, nday, nhour, 0 )
call Init_nmdays( current_date, JUMPOVER29FEB)
!*********initialize grid parameters*********
if(MANUAL_GRID)then
! define the grid parameter manually (explicitely)
if(MasterProc)write(*,*)'DEFINING GRID MANUALLY!'
! must set all parameters... see example in version from 20151019 (or before)
else
! NOT MANUAL GRID
!check first if grid is defined in a separate file:
filename='Grid_Def.nc'
inquire(file=filename,exist=Grid_Def_exist)
Grid_Def_exist=Grid_Def_exist.and.Use_Grid_Def
if(Grid_Def_exist)then
if(MasterProc)write(*,*)'Found Grid_Def! ',trim(filename)
else
if(MasterProc.and.Use_Grid_Def)&
write(*,*)'Did not found Grid_Def ',trim(filename)
filename=date2string(meteo,startdate,mode='YMDH')
end if
if(MasterProc)write(*,*)'reading domain sizes from ',trim(filename)
call GetFullDomainSize(filename,IIFULLDOM,JJFULLDOM,KMAX_MET,projection)
KMAX_MID=0!initialize
open(IO_TMP,file=Vertical_levelsFile,action="read",iostat=ios)
if(ios==0)then
! define own vertical coordinates
if(MasterProc)&
write(*,*)'Define vertical levels from ',trim(Vertical_levelsFile)
read(IO_TMP,*)KMAX_MID
if(MasterProc)write(*,*)KMAX_MID, 'vertical levels '
External_Levels_Def=.true.
! Must use eta coordinates
if(MasterProc.and..not.USE_EtaCOORDINATES)&
write(*,*)'WARNING: using hybrid levels even if not asked to! '
USE_EtaCOORDINATES=.true.
else
if(MasterProc)write(*,*)'WARNING: could not open '//trim(Vertical_levelsFile)
External_Levels_Def=.false.
close(IO_TMP)
KMAX_MID=KMAX_MET
end if
KMAX_BND=KMAX_MID+1
allocate(A_bnd(KMAX_BND),B_bnd(KMAX_BND))
allocate(A_mid(KMAX_MID),B_mid(KMAX_MID))
allocate(dA(KMAX_MID),dB(KMAX_MID),dEta_i(KMAX_MID))
allocate(sigma_bnd(KMAX_BND),sigma_mid(KMAX_MID))
allocate(Eta_bnd(KMAX_BND),Eta_mid(KMAX_MID))
allocate(i_local(IIFULLDOM))
allocate(j_local(JJFULLDOM))
! set RUNDOMAIN default values where not defined
if(RUNDOMAIN(1)<1)RUNDOMAIN(1)=1
if(RUNDOMAIN(2)<1 .or. RUNDOMAIN(2)>IIFULLDOM) RUNDOMAIN(2)=IIFULLDOM
if(RUNDOMAIN(3)<1)RUNDOMAIN(3)=1
if(RUNDOMAIN(4)<1 .or. RUNDOMAIN(4)>JJFULLDOM) RUNDOMAIN(4)=JJFULLDOM
if(MasterProc)then
55 format(A,I5,A,I5)
write(*,55) 'FULLDOMAIN has sizes ',IIFULLDOM,' X ',JJFULLDOM
write(IO_LOG,55)'FULLDOMAIN has sizes ',IIFULLDOM,' X ',JJFULLDOM
write(*,55) 'RUNDOMAIN x coordinates from ',RUNDOMAIN(1),' to ',RUNDOMAIN(2)
write(IO_LOG,55)'RUNDOMAIN x coordinates from ',RUNDOMAIN(1),' to ',RUNDOMAIN(2)
write(*,55) 'RUNDOMAIN y coordinates from ',RUNDOMAIN(3),' to ',RUNDOMAIN(4)
write(IO_LOG,55)'RUNDOMAIN y coordinates from ',RUNDOMAIN(3),' to ',RUNDOMAIN(4)
end if
call find_poles(filename,Pole_Singular)
MIN_GRIDS=5
if(NPROC==NPROC_MPI)then
! partition into subdomains
call parinit(MIN_GRIDS,Pole_Singular) ! subdomains sizes and position
else
! partition into largesubdomains and subdomains
call parinit_groups(MIN_GRIDS,Pole_Singular) ! subdomains sizes and position
end if
call Alloc_MetFields(LIMAX,LJMAX,KMAX_MID,KMAX_BND,NMET)
if(ME_CALC>=0)then
call Alloc_GridFields(LIMAX,LJMAX,KMAX_MID,KMAX_BND)
else
call Alloc_GridFields(largeLIMAX,largeLJMAX,KMAX_MID,KMAX_BND)
end if
if(ME_CALC>=0)then
call Getgridparams(LIMAX,LJMAX,filename,cyclicgrid)
! defines i_fdom,j_fdom,i_local,j_local,Cyclicgrid,North_pole,Poles
! GRIDWIDTH_M, glon, glat, xm_i,xm_j,xm2,xmd,xmdji,xm2ji,gl_stagg,gb_stagg
! P0,A_bnd_met,B_bnd_met,A_bnd,B_bnd,A_mid,B_mid,sigma_mid,sigma_bnd
! for Stereographic projection:
! ref_latitude,fi,xp,yp,AN
! for lon lat projection:
! no additional parameters
! for Rotated_Spherical projection:
! grid_north_pole_latitude,grid_north_pole_longitude,x1_rot,y1_rot,dx_rot,dx_roti
else
call Getgridparams(largeLIMAX,largeLJMAX,filename,cyclicgrid)
end if
if(ios==0)close(IO_TMP)
end if
end subroutine GridRead
subroutine GetFullDomainSize(filename,IIFULLDOM,JJFULLDOM,KMAX,projection)
! Get input grid sizes
implicit none
character (len = *), intent(in) ::filename
integer, intent(out):: IIFULLDOM,JJFULLDOM,KMAX
character (len = *), intent(out) ::projection
integer :: status,ncFileID,idimID,jdimID, kdimID,timeDimID
integer :: GIMAX_file,GJMAX_file,KMAX_file,wrf_proj_code
real :: wrf_POLE_LAT=0.0
character (len = 30) ::MAP_PROJ_CHAR
integer :: NTime_Read
real :: TimesInDays(1000)
if(ME_MPI==0)then
print *,'Defining grid properties from ',trim(filename)
! open an existing netcdf dataset
status = nf90_open(path=trim(filename),mode=nf90_nowrite,ncid=ncFileID)
if(status/=nf90_noerr) then
print *,'not found',trim(filename)
call StopAll("GridValues: File not found")
end if
projection=''
status = nf90_get_att(ncFileID,nf90_global,"projection",projection)
if(status/=nf90_noerr) then
! WRF projection format
call check(nf90_get_att(ncFileID,nf90_global,"MAP_PROJ",wrf_proj_code))
if(.not.USE_WRF_MET_NAMES .and. MasterProc)write(*,*)'Assuming WRF metdata'
USE_WRF_MET_NAMES = .true.
select case(wrf_proj_code)
case(6)
status = nf90_get_att(ncFileID,nf90_global,"POLE_LAT",wrf_POLE_LAT)
if(status==nf90_noerr) then
write(*,*)"POLE_LAT", wrf_POLE_LAT
if(abs(wrf_POLE_LAT-90.0)<0.001)then
projection='lon lat'
else
projection='Rotated_Spherical'
end if
else
write(*,*)"POLE_LAT not found"
projection='lon lat'
end if
case(2)
projection='Stereographic'
case(1)
projection='lambert'
call check(nf90_get_att(ncFileID,nf90_global,"MAP_PROJ_CHAR",MAP_PROJ_CHAR))
write(*,*)"wrf projection: "//trim(MAP_PROJ_CHAR)
case default
call CheckStop("Projection not recognized")
end select
end if
! put into emep standard
if(trim(projection)=='Polar Stereographic')projection='Stereographic'
if(trim(projection)=='Rotated_Spherical'.or.trim(projection)=='rotated_spherical'&
.or.trim(projection)=='rotated_pole'.or.trim(projection)=='rotated_latitude_longitude')then
projection='Rotated_Spherical'
end if
if(trim(projection)=='lambert_conformal_conic'.or.trim(projection)=='Lambert Conformal')projection='lambert'
write(*,*)'projection: ',trim(projection)
! get dimensions id
if(trim(projection)=='Stereographic') then
status = nf90_inq_dimid(ncid=ncFileID, name="i", dimID=idimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="west_east", dimID=idimID))
status = nf90_inq_dimid(ncid=ncFileID, name="j", dimID=jdimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="south_north", dimID=jdimID))
elseif(trim(projection)=='lambert') then
status = nf90_inq_dimid(ncid=ncFileID, name="x", dimID=idimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="west_east", dimID=idimID))
status = nf90_inq_dimid(ncid=ncFileID, name="y", dimID=jdimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="south_north", dimID=jdimID))
write(*,*)'x y dimensions'
elseif(trim(projection)==trim('lon lat')) then
status=nf90_inq_dimid(ncid=ncFileID, name="lon", dimID=idimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="west_east", dimID=idimID))
status=nf90_inq_dimid(ncid=ncFileID, name="lat", dimID=jdimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="south_north", dimID=jdimID))
else
! write(*,*)'GENERAL PROJECTION ',trim(projection)
status=nf90_inq_dimid(ncid=ncFileID, name="i", dimID = idimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="west_east", dimID=idimID))
status=nf90_inq_dimid(ncid=ncFileID, name="j", dimID = jdimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="south_north", dimID=jdimID))
end if
status=nf90_inq_dimid(ncid=ncFileID, name="k", dimID=kdimID)
if(status/=nf90_noerr)then
status=nf90_inq_dimid(ncid=ncFileID, name="lev", dimID=kdimID)!hybrid coordinates
if(status/=nf90_noerr) then
status=nf90_inq_dimid(ncid=ncFileID, name="hybrid", dimID=kdimID)!hybrid coordinates
if(status/=nf90_noerr) then ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="bottom_top", dimID=kdimID))
end if
end if
end if
!get dimensions length
call check(nf90_inquire_dimension(ncid=ncFileID,dimID=idimID,len=GIMAX_file))
call check(nf90_inquire_dimension(ncid=ncFileID,dimID=jdimID,len=GJMAX_file))
call check(nf90_inquire_dimension(ncid=ncFileID,dimID=kdimID,len=KMAX_file))
write(*,*)'dimensions input grid:',GIMAX_file,GJMAX_file,KMAX_file!,Nhh
IIFULLDOM=GIMAX_file
JJFULLDOM=GJMAX_file
KMAX =KMAX_file
call check(nf90_close(ncFileID))
end if
CALL MPI_BCAST(USE_WRF_MET_NAMES ,1,MPI_LOGICAL,0,MPI_COMM_WORLD,IERROR)
CALL MPI_BCAST(IIFULLDOM ,4*1,MPI_BYTE,0,MPI_COMM_WORLD,IERROR)
CALL MPI_BCAST(JJFULLDOM ,4*1,MPI_BYTE,0,MPI_COMM_WORLD,IERROR)
CALL MPI_BCAST(KMAX ,4*1,MPI_BYTE,0,MPI_COMM_WORLD,IERROR)
CALL MPI_BCAST(projection ,len(projection),MPI_BYTE,0,MPI_COMM_WORLD,IERROR)
end subroutine GetFullDomainSize
subroutine find_poles(filename,Pole_Singular)
! defines if there is a singularity at poles
implicit none
character (len = *), intent(in) ::filename
integer, intent(out):: Pole_Singular
integer :: status,ncFileID,varid
real,allocatable :: latitudes(:)
Pole_Singular=0
if(trim(projection)==trim('lon lat')) then
! find wether poles are included (or almost included) in grid
!
! If some cells are to narrow (Poles in lat lon coordinates),
! this will give too small time steps in the Advection,
! because of the constraint that the Courant number should be <1.
!
! If Poles are found and lon-lat coordinates are used the Advection scheme
! will be modified to be able to cope with the singularity
! the advection routine will not work efficiently with NPROCY>2 in this case
if(ME_MPI==0)then
! open an existing netcdf dataset
status = nf90_open(path=trim(filename),mode=nf90_nowrite,ncid=ncFileID)
if(status/=nf90_noerr) then
print *,'not found',trim(filename)
call StopAll("GridValues: File not found")
end if
allocate(latitudes(JJFULLDOM))
status=nf90_inq_varid(ncid=ncFileID, name="lat", varID=varID)
if(status/=nf90_noerr) then
!WRF format
call check(nf90_inq_varid(ncid=ncFileID, name="XLAT", varID=varID))
call check(nf90_get_var(ncFileID, varID,latitudes ,start=(/1,1/), count=(/1,JJFULLDOM/) ))
else
call check(nf90_get_var(ncFileID, varID,latitudes ))
endif
if(latitudes(RUNDOMAIN(4))>88.0)then
write(*,*)'The grid is singular at North Pole'
Pole_Singular=Pole_Singular+1
end if
if(latitudes(RUNDOMAIN(3))<-88.0)then
write(*,*)'The grid is singular at South Pole'
Pole_Singular=Pole_Singular+1
end if
deallocate(latitudes)
call check(nf90_close(ncFileID))
end if
end if
CALL MPI_BCAST(Pole_Singular ,4*1,MPI_BYTE,0,MPI_COMM_WORLD,IERROR)
end subroutine find_poles
subroutine Getgridparams(LIMAX,LJMAX,filename,cyclicgrid)
! Get grid and time parameters as defined in the meteo or Grid_Def file
! Do some checks on sizes and dates
!
! This routine is called only once (and is therefore not optimized for speed)
! defines i_fdom,j_fdom,i_local,j_local,Cyclicgrid,North_pole,Poles
! GRIDWIDTH_M, glon, glat, xm_i,xm_j,xm2,xmd,xmdji,xm2ji,gl_stagg,gb_stagg
! P0,A_bnd_met,B_bnd_met,A_bnd,B_bnd,A_mid,B_mid,sigma_mid,sigma_bnd
! for Stereographic projection:
! ref_latitude,fi,xp,yp,AN
! for lon lat projection:
! no additional parameters
! for Rotated_Spherical projection:
! grid_north_pole_latitude,grid_north_pole_longitude,x1_rot,y1_rot,dx_rot,dx_roti
implicit none
integer, intent(in):: LIMAX,LJMAX
character (len = *), intent(in) ::filename
integer, intent(out):: cyclicgrid
integer :: n,i,j,k,kk
integer :: ncFileID,idimID,jdimID,varID
integer :: status,South_pole,North_pole
real :: x1,x2,x3,x4,P0,x,y,mpi_out,r,t
logical::found_hybrid=.false.,found_metlevels=.false.
real :: CEN_LAT, CEN_LON,P_TOP_MET, WRF_DY
real :: rb,rl,rp,dx,dy,dy2,glmax,glmin,v2(2),glon_fdom1,glat_fdom1,lat
integer :: iloc_start, iloc_end,jloc_start, jloc_end
real, dimension(-1:LIMAX+2,-1:LJMAX+2)::xm,xm_i_ext,xm_j_ext
real, dimension(0:LIMAX+1,0:LJMAX+1)::lon_ext,lat_ext
!define longitudes in interval [-180,180]
glmin = -180.0
glmax = glmin + 360.0
! we can already define some arrays:
!/ Define full-domain coordinates of local i,j values. We need to account for
! the fact that each parallel domain has its starting cordinate
! gi0, gj0, and the user may specify a set of lower-left starting
! coordinates for running the model, IRUNBEG, JRUNBEG
! i_fdom(i) = i + gi0 + IRUNBEG - 2
! j_fdom(j) = j + gj0 + JRUNBEG - 2
i_fdom = (/ (n + gi0 + IRUNBEG - 2, n=0,LIMAX+1) /)
j_fdom = (/ (n + gj0 + JRUNBEG - 2, n=0,LJMAX+1) /)
! And the reverse, noting that we even define for area
! outside local domain
i_local = (/ (n - gi0 - IRUNBEG + 2, n=1, IIFULLDOM) /)
j_local = (/ (n - gj0 - JRUNBEG + 2, n=1, JJFULLDOM) /)
call CheckStop(GIMAX+IRUNBEG-1 > IIFULLDOM, "GridRead: I outside domain" )
call CheckStop(GJMAX+JRUNBEG-1 > JJFULLDOM, "GridRead: J outside domain" )
!open an existing netcdf dataset
status = nf90_open(path=trim(filename),mode=nf90_nowrite,ncid=ncFileID)
if(status/=nf90_noerr) then
print *,'not found',trim(filename)
call CheckStop("GridValues: File not found")
end if
if(MasterProc)print *,'Defining grid parameters from ',trim(filename)
!get dimensions id
select case(projection)
case('Stereographic')
status = nf90_inq_dimid(ncid=ncFileID, name="i", dimID=idimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="west_east", dimID=idimID))
status = nf90_inq_dimid(ncid=ncFileID, name="j", dimID=jdimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="south_north", dimID=jdimID))
case('lambert')
status = nf90_inq_dimid(ncid=ncFileID, name="x", dimID=idimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="west_east", dimID=idimID))
status = nf90_inq_dimid(ncid=ncFileID, name="y", dimID=jdimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="south_north", dimID=jdimID))
case('lon lat')
status=nf90_inq_dimid(ncid=ncFileID, name="lon", dimID=idimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="west_east", dimID=idimID))
status=nf90_inq_dimid(ncid=ncFileID, name="lat", dimID=jdimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="south_north", dimID=jdimID))
case default
! write(*,*)'GENERAL PROJECTION ',trim(projection)
status=nf90_inq_dimid(ncid=ncFileID, name="i", dimID=idimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="west_east", dimID=idimID))
status=nf90_inq_dimid(ncid=ncFileID, name="j", dimID=jdimID)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_inq_dimid(ncid=ncFileID, name="south_north", dimID=jdimID))
end select
!get global attributes
status = nf90_get_att(ncFileID,nf90_global,"Grid_resolution",GRIDWIDTH_M)
if(status/=nf90_noerr)then
!WRF format
call check(nf90_get_att(ncFileID,nf90_global,"DX",GRIDWIDTH_M))
status = nf90_get_att(ncFileID,nf90_global,"DY",v(1))
if(status==nf90_noerr .and. abs(GRIDWIDTH_M-v(1))>0.01) then
! if(MasterProc)write(*,*)'Gridcells not square. Will correct y mapping factor'
endif
end if
if(MasterProc)write(*,*)"Grid_resolution",GRIDWIDTH_M
select case(projection)
case('Stereographic')
status=nf90_get_att(ncFileID,nf90_global,"ref_latitude",ref_latitude)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_get_att(ncFileID,nf90_global,"TRUELAT1",ref_latitude))
status = nf90_get_att(ncFileID, nf90_global, "fi",fi )
if(status/=nf90_noerr)& ! WRF format
call check(nf90_get_att(ncFileID, nf90_global, "STAND_LON",fi ))
status = nf90_get_att(ncFileID, nf90_global, "xcoordinate_NorthPole",xp )
if(status==nf90_noerr) then
call check(nf90_get_att(ncFileID, nf90_global, "ycoordinate_NorthPole",yp ))
else
!WRF format compute from grid center coordinates
call check(nf90_get_att(ncFileID, nf90_global, "CEN_LAT", CEN_LAT))
call check(nf90_get_att(ncFileID, nf90_global, "CEN_LON", CEN_LON))
xp = 0.5 + 0.5*IIFULLDOM &
- EARTH_RADIUS/GRIDWIDTH_M*(1+sin(ref_latitude*PI/180.))&
*tan(PI/4-CEN_LAT*PI/180./2)*sin((CEN_LON-fi)*PI/180.)
yp = 0.5 + 0.5*JJFULLDOM &
+ EARTH_RADIUS/GRIDWIDTH_M*(1+sin(ref_latitude*PI/180.))&
*tan(PI/4-CEN_LAT*PI/180./2)*cos((CEN_LON-fi)*PI/180.)
!correct for last digits. Assume that numbers are close enough to an integer
if(abs(nint(xp)-xp)<0.01)xp=nint(xp)
if(abs(nint(yp)-yp)<0.01)yp=nint(yp)
if(MasterProc)then
write(*,*)"M= ",EARTH_RADIUS/GRIDWIDTH_M*(1+sin(ref_latitude*PI/180.))
write(*,*)"coordinates of North pole ",xp,yp
end if
end if
AN = 6.370e6*(1.0+sin(ref_latitude*PI/180.))/GRIDWIDTH_M
! = 237.7316364 for GRIDWIDTH_M=50 km and ref_latitude=60
do j = 0, LJMAX+1
dy = yp - j_fdom(j)
dy2 = dy*dy
do i = 0, LIMAX+1
dx = i_fdom(i) - xp
rp = sqrt(dx*dx+dy2) ! => distance to pole
rb = 90.0 - 180.0/PI*2* atan(rp/AN) ! => latitude
rl = 0.0
if (rp > 1.0e-10) rl = fi + 180.0/PI*atan2(dx,dy)
if (rl < glmin) rl = rl + 360.0
if (rl > glmax) rl = rl - 360.0
lon_ext(i,j)=rl ! longitude
lat_ext(i,j)=rb ! latitude
end do ! i
end do ! j
case('lambert')
if(USE_WRF_MET_NAMES)then
call check(nf90_get_att(ncFileID, nf90_global, "TRUELAT1",lat_stand1_lambert ))
call check(nf90_get_att(ncFileID, nf90_global, "TRUELAT2",lat_stand2_lambert ))
lat0_lambert = lat_stand1_lambert
call check(nf90_get_att(ncFileID, nf90_global, "STAND_LON",lon0_lambert ))
earth_radius_lambert = earth_radius
else
status = nf90_get_att(ncFileID,nf90_global,"latitude_of_projection_origin",lat0_lambert)!reference latitude, at which x=0
if(status/=nf90_noerr) then
call check(nf90_inq_varid(ncid=ncFileID, name="projection_lambert", varID=varID))
call check(nf90_get_att(ncFileID,varID,"latitude_of_projection_origin",lat0_lambert))
endif
status = nf90_get_att(ncFileID,nf90_global,"longitude_of_central_meridian",lon0_lambert)!reference longitude
if(status/=nf90_noerr) then
call check(nf90_inq_varid(ncid=ncFileID, name="projection_lambert", varID=varID))
call check(nf90_get_att(ncFileID,varID,"longitude_of_central_meridian",lon0_lambert))
endif
status = nf90_get_att(ncFileID,nf90_global,"earth_radius",earth_radius_lambert)!
if(status/=nf90_noerr) then
call check(nf90_inq_varid(ncid=ncFileID, name="projection_lambert", varID=varID))
call check(nf90_get_att(ncFileID,varID,"earth_radius",earth_radius_lambert))
endif
!status = nf90_get_att(ncFileID,nf90_global,"standard_parallel",(/lat_stand1_lambert,lat_stand2_lambert/))!standard latitude at which mapping factor=1
!default lat_stand1_lambert=lat_stand2_lambert=lat0_lambert
lat_stand1_lambert = lat0_lambert
lat_stand2_lambert = lat0_lambert
endif
if(abs(lat_stand1_lambert-lat_stand2_lambert)>1.E-4)then
k_lambert = log(cos(deg2rad*lat_stand1_lambert)/cos(deg2rad*lat_stand2_lambert))/&
(log(tan(0.25*PI+0.5*deg2rad*lat_stand2_lambert)/tan(0.25*PI+0.5*deg2rad*lat_stand1_lambert)))
lat0_lambert = rad2deg*asin(k_lambert)
if(MasterProc)then
write(*,*)'first true latitude ',lat_stand1_lambert
write(*,*)'second true latitude ',lat_stand2_lambert
write(*,*)'latitude of projection origin calculated to ',lat0_lambert
end if
else
k_lambert = sin(deg2rad*lat0_lambert)! also called n
endif
F_lambert = cos(deg2rad*lat_stand1_lambert) &
* tan(0.25*PI+0.5*deg2rad*lat_stand1_lambert)**k_lambert /k_lambert!normalization constant
y0_lambert = F_lambert*tan(0.25*PI-0.5*deg2rad*lat0_lambert)**k_lambert!reference y coordinate, also called rho0
if(USE_WRF_MET_NAMES)then
call check(nf90_inq_varid(ncid=ncFileID, name="XLONG", varID=varID))
call check(nf90_get_var(ncFileID,varID,v,count=(/1/)))
call check(nf90_inq_varid(ncid=ncFileID, name="XLAT", varID=varID))
call check(nf90_get_var(ncFileID,varID,u,count=(/1/)))
else
status = nf90_inq_varid(ncid=ncFileID, name="lon", varID=varID)
if(status/=nf90_noerr)&
call check(nf90_inq_varid(ncid=ncFileID, name="longitude", varID=varID))
call check(nf90_get_var(ncFileID,varID,v,count=(/1/)))
x1_lambert=v(1)
status = nf90_inq_varid(ncid=ncFileID, name="lat", varID=varID)
if(status/=nf90_noerr)&
call check(nf90_inq_varid(ncid=ncFileID, name="latitude", varID=varID))
call check(nf90_get_var(ncFileID,varID,u,count=(/1/)))
y1_lambert=u(1)
endif
x1_lambert=0.0
y1_lambert=0.0
call lb2ij(v(1),u(1),x1_lambert,y1_lambert)
x1_lambert = (x1_lambert-1)*GRIDWIDTH_M
y1_lambert = (y1_lambert-1)*GRIDWIDTH_M
if(MasterProc)then
! test that lon lat from meteo is the same as calculated, for the point (i,j)=(1,1)
call lb2ij(v(1),u(1),v(2),u(2))
if(abs(u(2)-1.0)+abs(v(2)-1.0)>1.E-4)then
write(*,*)'ERROR Lambert projection lon lat of (i,j)=(1,1) in file is',v(1),u(1)
write(*,*)'BUT Lambert projection (i,j) for this lon lat is not (1,1) but ',v(2),u(2)
call StopAll('ERROR in Lambert projection parameters')
else
write(*,*)'Lambert projection checked. (i,j)=(1,1) has lon lat',v(1),u(1)
write(*,*)"x and y at (i,j)=(1,1)",x1_lambert,y1_lambert
write(*,*)"y0_lambert,F_lambert ",y0_lambert,F_lambert
endif
endif
!make lon lat and mapping factors
do j = 0, LJMAX+1
y = (y1_lambert+(j_fdom(j)-1)*GRIDWIDTH_M)/EARTH_RADIUS
do i = 0, LIMAX+1
x = (x1_lambert+(i_fdom(i)-1)*GRIDWIDTH_M)/EARTH_RADIUS
r = sqrt(x*x+(y0_lambert-y)*(y0_lambert-y))
if(k_lambert<0.0)r = -r
t = atan(x/(y0_lambert-y))
lat_ext(i,j) = 2*rad2deg*atan((F_lambert/r)**(1.0/k_lambert))-90.0
lon_ext(i,j) = lon0_lambert + rad2deg*t/k_lambert
!does not work for lat = -90.0
xm(i,j)=k_lambert*F_lambert&
*tan(PI*0.25-deg2rad*0.5*lat_ext(i,j))**(k_lambert-1)&
*0.5/(cos(PI*0.25-deg2rad*0.5*lat_ext(i,j))**2)
xm2(i,j) = xm(i,j)*xm(i,j)
xmd(i,j) = 1.0/xm2(i,j)
xm2ji(j,i) = xm2(i,j)
xmdji(j,i) = xmd(i,j)
enddo
enddo
!staggered map factors
do j = 0, LJMAX+1
y = (y1_lambert+(j_fdom(j)-1+0.5)*GRIDWIDTH_M)/EARTH_RADIUS
do i = 0, LIMAX+1
x = (x1_lambert+(i_fdom(i)-1)*GRIDWIDTH_M)/EARTH_RADIUS
r = sqrt(x*x+(y0_lambert-y)*(y0_lambert-y))
if(k_lambert<0.0)r = -r
lat = 2*rad2deg*atan((F_lambert/r)**(1.0/k_lambert))-90.0
xm_i(i,j)=k_lambert*F_lambert&
*tan(PI*0.25-deg2rad*0.5*lat)**(k_lambert-1)&
*0.5/(cos(PI*0.25-deg2rad*0.5*lat)**2)
enddo
enddo
do j = 0, LJMAX+1
y = (y1_lambert+(j_fdom(j)-1)*GRIDWIDTH_M)/EARTH_RADIUS
do i = 0, LIMAX+1
x = (x1_lambert+(i_fdom(i)-1+0.5)*GRIDWIDTH_M)/EARTH_RADIUS
r = sqrt(x*x+(y0_lambert-y)*(y0_lambert-y))
if(k_lambert<0.0)r = -r
lat = 2*rad2deg*atan((F_lambert/r)**(1.0/k_lambert))-90.0
xm_j(i,j)=k_lambert*F_lambert&
*tan(PI*0.25-deg2rad*0.5*lat)**(k_lambert-1)&
*0.5/(cos(PI*0.25-deg2rad*0.5*lat)**2)
enddo
enddo
case('lon lat')
if(.not. USE_WRF_MET_NAMES)then
!NB: lon and lat are stored as 1 dimensional arrays
call check(nf90_inq_varid(ncid=ncFileID, name="lon", varID=varID))
call check(nf90_get_var(ncFileID, varID, lon_ext(1:limax,1),&
start=(/gi0+IRUNBEG-1/),count=(/limax/) ))
if(LIMAX>limax)&
lon_ext(LIMAX,1)=lon_ext(limax,1)+(lon_ext(limax,1)-lon_ext(limax-1,1))
lon_ext(0,1)=2*lon_ext(1,1)-lon_ext(2,1)
lon_ext(LIMAX+1,1)=2*lon_ext(LIMAX,1)-lon_ext(LIMAX-1,1)
do j=0,LJMAX+1
lon_ext(:,j)=lon_ext(:,1)
end do
call check(nf90_inq_varid(ncid=ncFileID,name="lat",varID=varID))
call check(nf90_get_var(ncFileID, varID, lat_ext(1,1:ljmax),&
start=(/gj0+JRUNBEG-1/),count=(/ljmax/) ))
lat_ext(1,LJMAX)=min(90.0,lat_ext(1,LJMAX))!should never be used anyway
lat_ext(1,0)=2*lat_ext(1,1)-lat_ext(1,2)
lat_ext(1,LJMAX+1)=2*lat_ext(1,LJMAX)-lat_ext(1,LJMAX-1)
do i=0,LIMAX+1
lat_ext(i,:)=lat_ext(1,:)
end do
else
!WRF format
call check(nf90_inq_varid(ncid=ncFileID, name="XLONG", varID=varID))
call nf90_get_var_extended(ncFileID,varID,lon_ext,0,LIMAX+1,0,LJMAX+1)
call check(nf90_inq_varid(ncid=ncFileID, name="XLAT", varID=varID))
call nf90_get_var_extended(ncFileID,varID,lat_ext,0,LIMAX+1,0,LJMAX+1)
end if
case('Rotated_Spherical')
status=nf90_get_att(ncFileID,nf90_global,"grid_north_pole_latitude",grid_north_pole_latitude)
if(status/=nf90_noerr)& ! WRF format
call check(nf90_get_att(ncFileID,nf90_global,"POLE_LAT",grid_north_pole_latitude))
if(MasterProc)write(*,*)"grid_north_pole_latitude",grid_north_pole_latitude
status=nf90_get_att(ncFileID,nf90_global,"grid_north_pole_longitude",grid_north_pole_longitude)
if(status/=nf90_noerr) then
! WRF format
call check(nf90_get_att(ncFileID,nf90_global,"POLE_LON",grid_north_pole_longitude))
!find resolution in degrees from resolution in km. WRF uses Erath Radius 6370 km(?)
dx_rot=360./(6370000.*2*PI/GRIDWIDTH_M)
!round to 6 digits
dx_rot=0.000001*nint(1000000*dx_rot)
end if
if(MasterProc)write(*,*)"grid_north_pole_longitude",grid_north_pole_longitude
status=nf90_inq_varid(ncid=ncFileID, name="i", varID=varID)
if(status==nf90_noerr) then
call check(nf90_get_var(ncFileID, varID, v2))!note that i is one dimensional
x1_rot=v2(1)
dx_rot=v2(2)-v2(1)
call check(nf90_inq_varid(ncid=ncFileID, name="j", varID=varID))
call check(nf90_get_var(ncFileID, varID, v2(1)))!note that j is one dimensional
y1_rot=v2(1)
call check(nf90_inq_varid(ncid=ncFileID, name="lon", varID=varID))
call nf90_get_var_extended(ncFileID,varID,lon_ext,0,LIMAX+1,0,LJMAX+1)
call check(nf90_inq_varid(ncid=ncFileID, name="lat", varID=varID))
call nf90_get_var_extended(ncFileID,varID,lat_ext,0,LIMAX+1,0,LJMAX+1)
else
! WRF format
call check(nf90_inq_varid(ncid=ncFileID, name="XLONG", varID=varID))
call nf90_get_var_extended(ncFileID,varID,lon_ext,0,LIMAX+1,0,LJMAX+1)
call check(nf90_get_var(ncFileID, varID, v2,start=(/1,1/),count=(/1,1/) ))
glon_fdom1=v2(1)
! glon=0.0!to get some value for outside subdomain too (when limax<LIMAX for instance)
! call check(nf90_get_var(ncFileID, varID, glon(1:limax,1:ljmax),&
! start=(/gi0+IRUNBEG-1,gj0+JRUNBEG-1/),count=(/limax,ljmax/) ))
call check(nf90_inq_varid(ncid=ncFileID, name="XLAT", varID=varID))
call nf90_get_var_extended(ncFileID,varID,lat_ext,0,LIMAX+1,0,LJMAX+1)
call check(nf90_get_var(ncFileID, varID, v2,start=(/1,1/),count=(/1,1/) ))
glat_fdom1=v2(1)
x1_rot=0.
y1_rot=0.
dx_roti=1.0/dx_rot
call lb2ij(glon_fdom1,glat_fdom1,x,y)
x1_rot=(x-1)*dx_rot
y1_rot=(y-1)*dx_rot
do i=1,10
if(x1_rot>180.0)then
x1_rot=x1_rot-360.0
else if(x1_rot<-180.0)then
x1_rot=x1_rot+360.0
else
exit
end if
end do
! call lb2ij(glon_fdom(1,1),glat_fdom(1,1),x,y)
! write(*,*)'after ',glon_fdom(1,1),glat_fdom(1,1),x,y
! call lb_rot2lb(x,y,x1_rot,y1_rot,grid_north_pole_longitude,grid_north_pole_latitude)
! write(*,*)"spherical lon lat of (i,j)=(1,1)",x,y,glon_fdom(1,1),glat_fdom(1,1)
if(MasterProc)write(*,*)"rotated lon lat of (i,j)=(1,1)",x1_rot,y1_rot
if(MasterProc)write(*,*)"resolution",dx_rot
end if
dx_roti=1.0/dx_rot
case default
! other projection?
call check(nf90_inq_varid(ncid=ncFileID, name="lon", varID=varID))
call nf90_get_var_extended(ncFileID,varID,lon_ext,0,LIMAX+1,0,LJMAX+1)
call check(nf90_inq_varid(ncid=ncFileID, name="lat", varID=varID))
call nf90_get_var_extended(ncFileID,varID,lat_ext,0,LIMAX+1,0,LJMAX+1)
end select
glon(1:LIMAX,1:LJMAX)=lon_ext(1:LIMAX,1:LJMAX) ! longitude
glat(1:LIMAX,1:LJMAX)=lat_ext(1:LIMAX,1:LJMAX) ! latitude
do j=1,LJMAX
do i=1,LIMAX
if(glon(i,j)>glmax)glon(i,j)=glon(i,j)-360.0
if(glon(i,j)<glmin)glon(i,j)=glon(i,j)+360.0
end do
end do
if(projection=='lon lat')then
!we require the longitudes to increase if i increases
if(glon(1,1)>glon(limax,1))then
write(*,*)'WARNING: shifting longitudes by 360 deg for processor ',me,&
' lon was from ', glon(1,1),'to ',glon(limax,1)
do j=1,LJMAX
do i=1,LIMAX
if(glon(i,j)<glon(1,j))glon(i,j)=glon(i,j)+360.0
end do
end do
endif
endif
! map factors
status=nf90_inq_varid(ncid=ncFileID, name="map_factor", varID=varID)
if(status==nf90_noerr)then
! mapping factor at center of cells is defined
call nf90_get_var_extended(ncFileID,varID,xm,-1,LIMAX+2,-1,LJMAX+2)
!make "staggered" map factors and other derived map factors
do j=0,LJMAX+1
do i=0,LIMAX+1
xm_i(i,j)=0.5*(xm(i,j)+xm(i,j+1))
xm_j(i,j)=0.5*(xm(i,j)+xm(i+1,j))
xm2(i,j)=xm(i,j)*xm(i,j)
xmd(i,j) =1.0/xm2(i,j)
xm2ji(j,i) = xm2(i,j)
xmdji(j,i) = xmd(i,j)
end do
end do
elseif(trim(projection)=='lambert') then
! map factors have been computed already (above)
else
!map factor are already staggered
status=nf90_inq_varid(ncid=ncFileID, name="map_factor_i", varID=varID)
iloc_start=-1
if(iloc_start+IRUNBEG+gi0-2<1)iloc_start=1!first cell (in i direction)
iloc_end=LIMAX+2
if(iloc_end+IRUNBEG+gi0-2>IIFULLDOM)iloc_end=IIFULLDOM+2-gi0-IRUNBEG!last cell
jloc_start=-1
if(jloc_start+JRUNBEG+gj0-2<1)jloc_start=1!first cell (in j direction)
jloc_end=LJMAX+2
if(jloc_end+JRUNBEG+gj0-2>JJFULLDOM)jloc_end=JJFULLDOM+2-gj0-JRUNBEG!last cell
if(status==nf90_noerr)then
call nf90_get_var_extended(ncFileID,varID,xm_i_ext,-1,LIMAX+2,-1,LJMAX+2)
else
!WRF format
call check(nf90_inq_varid(ncid=ncFileID, name="MAPFAC_VX", varID=varID))
call nf90_get_var_extended(ncFileID,varID,xm_i_ext,-1,LIMAX+2,-1,LJMAX+2,&
jshift_in=1) !NB:shift j by 1 since wrf start at bottom face
end if
status=nf90_inq_varid(ncid=ncFileID, name="map_factor_j", varID=varID)
if(status==nf90_noerr)then
call nf90_get_var_extended(ncFileID,varID,xm_j_ext,-1,LIMAX+2,-1,LJMAX+2)
else
!WRF format
call check(nf90_inq_varid(ncid=ncFileID, name="MAPFAC_UY", varID=varID))
call nf90_get_var_extended(ncFileID,varID,xm_j_ext,-1,LIMAX+2,-1,LJMAX+2,&
ishift_in=1) !NB:shift i by 1 since wrf start at left face
status = nf90_get_att(ncFileID,nf90_global,"DY",WRF_DY)
if(status==nf90_noerr) then
!WRF uses DY/MAPFAC_UY while emep uses GRIDWIDTH_M/xm_j for the y size
if(abs(GRIDWIDTH_M/WRF_DY-1.0)>1.E-6)then
if(MasterProc)write(*,*)"rescaling y mapfactors with = ",GRIDWIDTH_M/WRF_DY
xm_j_ext=xm_j_ext*GRIDWIDTH_M/WRF_DY
endif
else
if(MasterProc)write(*,*)"not rescaling y mapfactors"
endif
end if
!define xm2, xm_i and xm_j now
!Note that xm is inverse length: interpolate 1/xm rather than xm
do j=0,LJMAX+1
do i=0,LIMAX+1
xm_i(i,j)=max(1.0E-5,xm_i_ext(i,j))
xm_j(i,j)=max(1.0E-5,xm_j_ext(i,j))
xm2(i,j) = 4.0*( (xm_i_ext(i,j-1)*xm_i_ext(i,j))/&
(xm_i_ext(i,j-1)+xm_i_ext(i,j)) )&
*( (xm_j_ext(i-1,j)*xm_j_ext(i,j))/&
(xm_j_ext(i-1,j)+xm_j_ext(i,j)) )
xm2(i,j)=max(1.E-7,xm2(i,j))
xmd(i,j) =1.0/xm2(i,j)
xm2ji(j,i) = xm2(i,j)
xmdji(j,i) = xmd(i,j)
end do
end do
end if
status=nf90_inq_varid(ncid=ncFileID, name="k", varID=varID)