-
Notifications
You must be signed in to change notification settings - Fork 321
/
Copy pathlilac_atmaero.F90
348 lines (291 loc) · 14.6 KB
/
lilac_atmaero.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
module lilac_atmaero
!-----------------------------------------------------------------------
! Contains methods for reading in atmosphere aerosal data
! This will be done on the CTSM grid with the CTSM decomposition
! (after the redistribution from atm-> lnd)
!-----------------------------------------------------------------------
use ESMF
! share code uses
use shr_kind_mod , only : r8 => shr_kind_r8, CL => shr_kind_cl, CS => shr_kind_cs
use shr_sys_mod , only : shr_sys_abort
use shr_nl_mod , only : shr_nl_find_group_name
use shr_log_mod , only : shr_log_errMsg
use shr_mpi_mod , only : shr_mpi_bcast
use shr_strdata_mod , only : shr_strdata_type, shr_strdata_create
use shr_strdata_mod , only : shr_strdata_print, shr_strdata_advance
use shr_string_mod , only : shr_string_listAppend
use shr_cal_mod , only : shr_cal_ymd2date
use shr_pio_mod , only : shr_pio_getiotype
use mct_mod , only : mct_avect_indexra, mct_gsmap, mct_ggrid
use mct_mod , only : mct_gsmap_init, mct_gsmap_orderedpoints
use mct_mod , only : mct_ggrid_init, mct_ggrid_importIAttr, mct_ggrid_importRattr
! ctsm uses
use ncdio_pio , only : pio_subsystem
use ctsm_Domain , only : ldomain
use ctsm_TimeManager , only : get_calendar
! lilac uses
use lilac_atmcap , only : gindex_atm
use lilac_methods , only : chkerr
use lilac_methods , only : lilac_methods_FB_getFieldN
use lilac_constants , only : field_index_unset
use ctsm_LilacCouplingFields, only : a2l_fields, lilac_atm2lnd
use ctsm_LilacCouplingFieldIndices
implicit none
private
type, private :: field_mapping_type
character(len=:), allocatable :: field_name
integer :: field_index = field_index_unset
end type field_mapping_type
public :: lilac_atmaero_init ! initialize stream data type sdat
public :: lilac_atmaero_interp ! interpolates between two years of ndep file data
! module data
type(shr_strdata_type) :: sdat ! input data stream
! The first num_fields_to_read in the fields_to_read list are the fields that this
! module will read from data. This is set up to have the same ordering as the fields in
! sdat.
integer :: num_fields_to_read
type(field_mapping_type), allocatable :: fields_to_read(:)
character(*),parameter :: u_file_u = &
__FILE__
!==============================================================================
contains
!==============================================================================
subroutine lilac_atmaero_init(atm2cpl_state, rc)
! ----------------------------------------
! Initialize data stream information.
! ----------------------------------------
! input/output variables
type(ESMF_State) , intent(inout) :: atm2cpl_state
integer , intent(out) :: rc
! local variables
type(ESMF_VM) :: vm
type(ESMF_Mesh) :: lmesh
type(ESMF_FieldBundle) :: lfieldbundle
type(ESMF_Field) :: lfield
type(mct_ggrid) :: ggrid_atm ! domain information
type(mct_gsmap) :: gsmap_atm ! decompositoin info
type(field_mapping_type), allocatable :: all_fields(:) ! all fields that can possibly be read from data
integer :: mytask ! mpi task number
integer :: mpicom ! mpi communicator
integer :: n ! index
integer :: field_index
integer :: lsize ! local size
integer :: gsize ! global size
integer :: nunit ! namelist input unit
integer :: ierr ! namelist i/o error flag
character(len=cl) :: stream_fldfilename ! name of input stream file
character(len=CL) :: mapalgo = 'bilinear' ! type of 2d mapping
character(len=CS) :: taxmode = 'extend' ! time extrapolation
character(len=CL) :: fldlistFile ! name of fields in input stream file
character(len=CL) :: fldlistModel ! name of fields in model
integer :: stream_year_first ! first year in stream to use
integer :: stream_year_last ! last year in stream to use
integer :: model_year_align ! align stream_year_first with model year
integer :: spatialDim
integer :: numOwnedElements
real(r8), pointer :: ownedElemCoords(:)
real(r8), pointer :: mesh_lons(:)
real(r8), pointer :: mesh_lats(:)
real(r8), pointer :: mesh_areas(:)
real(r8), pointer :: rdata(:)
integer , pointer :: idata(:)
!-----------------------------------------------------------------------
namelist /atmaero_stream/ &
stream_year_first, stream_year_last, model_year_align, &
stream_fldfilename
rc = ESMF_SUCCESS
all_fields = [ &
field_mapping_type('BCDEPWET', lilac_a2l_Faxa_bcphiwet), &
field_mapping_type('BCPHODRY', lilac_a2l_Faxa_bcphodry), &
field_mapping_type('BCPHIDRY', lilac_a2l_Faxa_bcphidry), &
field_mapping_type('OCDEPWET', lilac_a2l_Faxa_ocphiwet), &
field_mapping_type('OCPHIDRY', lilac_a2l_Faxa_ocphidry), &
field_mapping_type('OCPHODRY', lilac_a2l_Faxa_ocphodry), &
field_mapping_type('DSTX01WD', lilac_a2l_Faxa_dstwet1), &
field_mapping_type('DSTX01DD', lilac_a2l_Faxa_dstdry1), &
field_mapping_type('DSTX02WD', lilac_a2l_Faxa_dstwet2), &
field_mapping_type('DSTX02DD', lilac_a2l_Faxa_dstdry2), &
field_mapping_type('DSTX03WD', lilac_a2l_Faxa_dstwet3), &
field_mapping_type('DSTX03DD', lilac_a2l_Faxa_dstdry3), &
field_mapping_type('DSTX04WD', lilac_a2l_Faxa_dstwet4), &
field_mapping_type('DSTX04DD', lilac_a2l_Faxa_dstdry4)]
num_fields_to_read = 0
allocate(fields_to_read(size(all_fields)))
fldlistFile = ' '
fldlistModel = ' '
do n = 1, size(all_fields)
field_index = all_fields(n)%field_index
if (a2l_fields%is_needed_from_data(field_index)) then
num_fields_to_read = num_fields_to_read + 1
fields_to_read(num_fields_to_read) = all_fields(n)
call shr_string_listAppend(fldlistFile, fields_to_read(num_fields_to_read)%field_name)
call shr_string_listAppend(fldlistModel, a2l_fields%get_fieldname(field_index))
end if
end do
if (num_fields_to_read == 0) then
return
end if
! default values for namelist
stream_year_first = 1 ! first year in stream to use
stream_year_last = 1 ! last year in stream to use
model_year_align = 1 ! align stream_year_first with this model year
stream_fldFileName = ' '
! get mytask and mpicom
call ESMF_VMGetCurrent(vm, rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
call ESMF_VMGet(vm, localPet=mytask, mpiCommunicator=mpicom, rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
! Read namelist
if (mytask == 0) then
open(newunit=nunit, file='lilac_in', status='old', iostat=ierr )
call shr_nl_find_group_name(nunit, 'atmaero_stream', status=ierr)
if (ierr == 0) then
read(nunit, atmaero_stream, iostat=ierr)
if (ierr /= 0) then
call shr_sys_abort(' ERROR reading namelist '//shr_log_errMsg(u_file_u, __LINE__))
end if
else
call shr_sys_abort(' ERROR finding namelist '//shr_log_errMsg(u_file_u, __LINE__))
end if
close(nunit)
endif
call shr_mpi_bcast(stream_year_first , mpicom)
call shr_mpi_bcast(stream_year_last , mpicom)
call shr_mpi_bcast(model_year_align , mpicom)
call shr_mpi_bcast(stream_fldfilename, mpicom)
if (mytask == 0) then
print *, ' '
print *, 'atmaero stream settings:'
print *, ' stream_year_first = ',stream_year_first
print *, ' stream_year_last = ',stream_year_last
print *, ' model_year_align = ',model_year_align
print *, ' stream_fldFileName = ',stream_fldFileName
print *, ' '
endif
! ------------------------------
! create the mct gsmap
! ------------------------------
lsize = size(gindex_atm)
gsize = ldomain%ni * ldomain%nj
call mct_gsmap_init( gsmap_atm, gindex_atm, mpicom, 1, lsize, gsize )
! ------------------------------
! obtain mesh lats, lons and areas
! ------------------------------
call ESMF_StateGet(atm2cpl_state, 'a2c_fb', lfieldbundle, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return
call lilac_methods_FB_getFieldN(lfieldbundle, fieldnum=1, field=lfield, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
call ESMF_FieldGet(lfield, mesh=lmesh, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
call ESMF_MeshGet(lmesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return
if (numOwnedElements /= lsize) then
call shr_sys_abort('ERROR: numOwnedElements is not equal to lsize')
end if
allocate(ownedElemCoords(spatialDim*numOwnedElements))
call ESMF_MeshGet(lmesh, ownedElemCoords=ownedElemCoords, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return
allocate(mesh_lons(numOwnedElements))
allocate(mesh_lats(numOwnedElements))
allocate(mesh_areas(numOwnedElements))
do n = 1,numOwnedElements
mesh_lons(n) = ownedElemCoords(2*n-1)
mesh_lats(n) = ownedElemCoords(2*n)
mesh_areas(n) = 1.e36 ! hard-wire for now for testing
end do
! ------------------------------
! create the mct ggrid
! ------------------------------
call mct_ggrid_init( ggrid=ggrid_atm, CoordChars='lat:lon:hgt', OtherChars='area:aream:mask:frac', lsize=lsize)
call mct_gsmap_orderedpoints(gsmap_atm, mytask, idata)
call mct_gGrid_importIAttr(ggrid_atm,'GlobGridNum', idata, lsize)
call mct_gGrid_importRattr(ggrid_atm,"lon" , mesh_lons , lsize)
call mct_gGrid_importRattr(ggrid_atm,"lat" , mesh_lats , lsize)
call mct_gGrid_importRattr(ggrid_atm,"area", mesh_areas, lsize)
allocate(rdata(lsize))
rdata(:) = 1._R8
call mct_gGrid_importRattr(ggrid_atm,"mask", rdata, lsize)
deallocate(mesh_lons, mesh_lats, mesh_areas, rdata)
! ------------------------------
! create the stream data sdat
! ------------------------------
call shr_strdata_create(sdat,&
name = "atmaero", &
pio_subsystem = pio_subsystem, &
pio_iotype = shr_pio_getiotype(compid= 1), &
mpicom = mpicom, &
compid = 1, &
gsmap = gsmap_atm, &
ggrid = ggrid_atm, &
nxg = ldomain%ni, &
nyg = ldomain%nj, &
yearFirst = stream_year_first, &
yearLast = stream_year_last, &
yearAlign = model_year_align, &
offset = 0, &
domFilePath = '', &
domfilename = trim(stream_fldfilename), &
domTvarName = 'time', &
domXvarName = 'lon' , &
domYvarName = 'lat' , &
domAreaName = 'area', &
domMaskName = 'mask', &
filePath = '', &
filename = (/trim(stream_fldfilename)/), &
fldListFile = trim(fldlistFile), &
fldListModel = trim(fldlistModel), &
fillalgo = 'none', &
mapalgo = mapalgo, &
calendar = get_calendar(), &
taxmode = taxmode )
if (mytask == 0) then
call shr_strdata_print(sdat,'ATMAERO data')
endif
end subroutine lilac_atmaero_init
!================================================================
subroutine lilac_atmaero_interp(clock, rc)
! input/output variables
type(ESMF_Clock) :: clock
integer, intent(out) :: rc
! local variables
type(ESMF_VM) :: vm
integer :: mpicom ! mpi communicator
integer :: mytask ! mpi task number
type(ESMF_FieldBundle) :: lfieldbundle
type(ESMF_Time) :: currTime
integer :: yy, mm, dd, sec, curr_ymd
integer :: n
character(len=*), parameter :: subname='lilac_atmaero: [lilac_atmaero_interp]'
!-----------------------------------------------------------------------
rc = ESMF_SUCCESS
if (num_fields_to_read == 0) then
return
end if
! get mytask and mpicom
call ESMF_VMGetCurrent(vm, rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
call ESMF_VMGet(vm, localPet=mytask, mpiCommunicator=mpicom, rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
! get current time info
call ESMF_ClockGet( clock, currTime=currTime, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return
call ESMF_TimeGet( currTime, yy=yy, mm=mm, dd=dd, s=sec, rc=rc )
if (ChkErr(rc,__LINE__,u_FILE_u)) return
call shr_cal_ymd2date(yy,mm,dd,curr_ymd)
! advance the streams
call shr_strdata_advance(sdat, curr_ymd, sec, mpicom, 'atmaero')
do n = 1, num_fields_to_read
call set_field(n)
end do
end subroutine lilac_atmaero_interp
!==============================================================================
subroutine set_field(fieldnum)
! input/output data
integer, intent(in) :: fieldnum ! index into fields_to_read and sdat (which are assumed to have the same ordering)
! local data
integer :: field_index ! index in a2l_fields
!-----------------------------------------------------------------------
field_index = fields_to_read(fieldnum)%field_index
call lilac_atm2lnd(field_index, sdat%avs(1)%rAttr(fieldnum,:))
end subroutine set_field
end module lilac_atmaero