Line data Source code
1 : !-------------------------------------------------------------------------------
2 : ! Initializes the CAM physics grid mesh
3 : !-------------------------------------------------------------------------------
4 : module edyn_phys_grid
5 : use shr_kind_mod, only: r8 => shr_kind_r8, cs=>shr_kind_cs, cl=>shr_kind_cl
6 : use cam_logfile, only: iulog
7 : use cam_abortutils, only: endrun
8 :
9 : implicit none
10 :
11 : private
12 :
13 : public :: edyn_phys_grid_init
14 :
15 : contains
16 :
17 768 : subroutine edyn_phys_grid_init()
18 : use ESMF, only: ESMF_DistGrid, ESMF_DistGridCreate, ESMF_MeshCreate
19 : use ESMF, only: ESMF_FILEFORMAT_ESMFMESH,ESMF_MeshGet,ESMF_Mesh
20 : use phys_control, only: phys_getopts
21 : use phys_grid, only: get_ncols_p, get_gcol_p, get_rlon_all_p, get_rlat_all_p
22 : use ppgrid, only: begchunk, endchunk
23 : use edyn_esmf, only: edyn_esmf_chkerr, edyn_esmf_update_phys_mesh
24 : use shr_const_mod,only: shr_const_pi
25 : use ppgrid, only: pcols
26 : use error_messages,only: alloc_err
27 :
28 : ! Local variables
29 : integer :: ncols
30 : integer :: chnk, col, dindex
31 768 : integer, allocatable :: decomp(:)
32 : character(len=cl) :: grid_file
33 : character(len=*), parameter :: subname = 'edyn_gcomp_init'
34 : real(r8) , parameter :: radtodeg = 180.0_r8/shr_const_pi
35 : integer :: spatialDim
36 : integer :: numOwnedElements
37 768 : real(r8), pointer :: ownedElemCoords(:)
38 768 : real(r8), pointer :: lat(:), latMesh(:)
39 768 : real(r8), pointer :: lon(:), lonMesh(:)
40 : real(r8) :: lats(pcols) ! array of chunk latitudes
41 : real(r8) :: lons(pcols) ! array of chunk longitude
42 : integer :: i, c, n
43 : character(len=cs) :: tempc1,tempc2
44 : character(len=300) :: errstr
45 :
46 : ! dist_grid_2d: DistGrid for 2D fields
47 : type(ESMF_DistGrid) :: dist_grid_2d
48 :
49 : ! phys_mesh: Local copy of physics grid
50 : type(ESMF_Mesh) :: phys_mesh
51 :
52 : real(r8), parameter :: abstol = 1.e-6_r8
53 : integer :: total_cols, rc
54 :
55 : ! Find the physics grid file
56 768 : call phys_getopts(physics_grid_out=grid_file)
57 : ! Compute the local decomp
58 768 : total_cols = 0
59 3072 : do chnk = begchunk, endchunk
60 3072 : total_cols = total_cols + get_ncols_p(chnk)
61 : end do
62 2304 : allocate(decomp(total_cols), stat=rc)
63 768 : call alloc_err(rc,subname,'decomp',total_cols)
64 :
65 768 : dindex = 0
66 3072 : do chnk = begchunk, endchunk
67 2304 : ncols = get_ncols_p(chnk)
68 30720 : do col = 1, ncols
69 27648 : dindex = dindex + 1
70 29952 : decomp(dindex) = get_gcol_p(chnk, col)
71 : end do
72 : end do
73 :
74 : ! Create a DistGrid based on the physics decomp
75 768 : dist_grid_2d = ESMF_DistGridCreate(arbSeqIndexList=decomp, rc=rc)
76 768 : call edyn_esmf_chkerr(subname, 'ESMF_DistGridCreate phys decomp', rc)
77 :
78 : ! Create an ESMF_mesh for the physics decomposition
79 : phys_mesh = ESMF_MeshCreate(trim(grid_file), ESMF_FILEFORMAT_ESMFMESH, &
80 768 : elementDistgrid=dist_grid_2d, rc=rc)
81 768 : call edyn_esmf_chkerr(subname, 'ESMF_MeshCreateFromFile', rc)
82 :
83 768 : call edyn_esmf_update_phys_mesh(phys_mesh)
84 :
85 : ! Check that the mesh coordinates are consistent with the model physics column coordinates
86 :
87 : ! obtain mesh lats and lons
88 768 : call ESMF_MeshGet(phys_mesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc)
89 768 : call edyn_esmf_chkerr(subname, 'ESMF_MeshGet', rc)
90 :
91 768 : if (numOwnedElements /= total_cols) then
92 0 : write(tempc1,'(i10)') numOwnedElements
93 0 : write(tempc2,'(i10)') total_cols
94 : call endrun(trim(subname)//": ERROR numOwnedElements "// &
95 0 : trim(tempc1) //" not equal to local size "// trim(tempc2))
96 : end if
97 :
98 2304 : allocate(ownedElemCoords(spatialDim*numOwnedElements), stat=rc)
99 768 : call alloc_err(rc,subname,'ownedElemCoords',spatialDim*numOwnedElements)
100 :
101 2304 : allocate(lonMesh(total_cols), stat=rc)
102 768 : call alloc_err(rc,subname,'lonMesh',total_cols)
103 :
104 2304 : allocate(latMesh(total_cols), stat=rc)
105 768 : call alloc_err(rc,subname,'latMesh',total_cols)
106 :
107 768 : call ESMF_MeshGet(phys_mesh, ownedElemCoords=ownedElemCoords)
108 :
109 28416 : do n = 1,total_cols
110 27648 : lonMesh(n) = ownedElemCoords(2*n-1)
111 28416 : latMesh(n) = ownedElemCoords(2*n)
112 : end do
113 :
114 : ! obtain internally generated cam lats and lons
115 2304 : allocate(lon(total_cols), stat=rc);
116 768 : call alloc_err(rc,subname,'lon',total_cols)
117 :
118 28416 : lon(:) = 0._r8
119 :
120 2304 : allocate(lat(total_cols), stat=rc);
121 768 : call alloc_err(rc,subname,'lat',total_cols)
122 :
123 28416 : lat(:) = 0._r8
124 :
125 768 : n=0
126 3072 : do c = begchunk, endchunk
127 2304 : ncols = get_ncols_p(c)
128 : ! latitudes and longitudes returned in radians
129 2304 : call get_rlat_all_p(c, ncols, lats)
130 2304 : call get_rlon_all_p(c, ncols, lons)
131 30720 : do i=1,ncols
132 27648 : n = n+1
133 27648 : lat(n) = lats(i)*radtodeg
134 29952 : lon(n) = lons(i)*radtodeg
135 : end do
136 : end do
137 :
138 768 : errstr = ''
139 : ! error check differences between internally generated lons and those read in
140 28416 : do n = 1,total_cols
141 27648 : if (abs(lonMesh(n) - lon(n)) > abstol) then
142 192 : if ( (abs(lonMesh(n)-lon(n)) > 360._r8+abstol) .or. (abs(lonMesh(n)-lon(n)) < 360._r8-abstol) ) then
143 0 : write(errstr,100) n,lon(n),lonMesh(n), abs(lonMesh(n)-lon(n))
144 0 : write(iulog,*) trim(errstr)
145 : endif
146 : end if
147 28416 : if (abs(latMesh(n) - lat(n)) > abstol) then
148 : ! poles in the 4x5 SCRIP file seem to be off by 1 degree
149 0 : if (.not.( (abs(lat(n))>88.0_r8) .and. (abs(latMesh(n))>88.0_r8) )) then
150 0 : write(errstr,101) n,lat(n),latMesh(n), abs(latMesh(n)-lat(n))
151 0 : write(iulog,*) trim(errstr)
152 : endif
153 : end if
154 : end do
155 :
156 768 : if ( len_trim(errstr) > 0 ) then
157 0 : call endrun(subname//': physics mesh coords do not match model coords')
158 : end if
159 :
160 : ! deallocate memory
161 768 : deallocate(ownedElemCoords)
162 768 : deallocate(lon, lonMesh)
163 768 : deallocate(lat, latMesh)
164 768 : deallocate(decomp)
165 :
166 : 100 format('edyn_gcomp_init: coord mismatch... n, lon(n), lonmesh(n), diff_lon = ',i6,2(f21.13,3x),d21.5)
167 : 101 format('edyn_gcomp_init: coord mismatch... n, lat(n), latmesh(n), diff_lat = ',i6,2(f21.13,3x),d21.5)
168 :
169 2304 : end subroutine edyn_phys_grid_init
170 :
171 :
172 : end module edyn_phys_grid
|
