FV3-JEDI
fv3jedi_geos_mod.F90
Go to the documentation of this file.
1 ! (C) Copyright 2017-2018 UCAR
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 
6 #define I_AM_MAIN
7 #include "MAPL_Generic.h"
8 #include "unused_dummy.H"
9 
10 module fv3jedi_geos_mod
11 
12 use iso_c_binding
13 use datetime_mod
14 use duration_mod
15 use netcdf
16 
17 use fckit_configuration_module, only: fckit_configuration
18 
19 use fv3jedi_geom_mod, only: fv3jedi_geom
20 use fv3jedi_field_mod, only: fv3jedi_field, field_clen
21 use fv3jedi_increment_mod, only: fv3jedi_increment
22 use fv3jedi_kinds_mod, only: kind_real
23 use fv3jedi_state_mod, only: fv3jedi_state
24 use fckit_mpi_module, only: fckit_mpi_comm
25 
26 use mpi
27 use esmf
28 use mapl_capmod, only: mapl_cap
29 use mapl_applicationsupport
30 
31 use mapl_mod
32 use mapl_profiler, only: get_global_time_profiler, baseprofiler, timeprofiler
33 use geos_gcmgridcompmod, only: geos_gcm => setservices
34 
35 implicit none
36 private
37 
38 public :: geos_model
39 public :: geos_create
40 public :: geos_delete
41 public :: geos_initialize
42 public :: geos_step
43 public :: geos_finalize
44 
45 ! ------------------------------------------------------------------------------
46 
47 !> Fortran derived type to hold model definition
48 type :: geos_model
49  type(mapl_cap) :: cap
50  integer :: geossubsteps
51  type(fckit_mpi_comm) :: f_comm
52  integer :: isc, iec, jsc, jec, npz
53  type(esmf_time) :: start_time
54  logical :: saved_state
55  logical :: reforecast
56 end type geos_model
57 
58 ! --------------------------------------------------------------------------------------------------
59 
60  type :: field_attributes
61  type(esmf_field) :: field(1)
62  character(len=ESMF_MAXSTR) :: short_name, long_name, units
63  end type field_attributes
64 
65 contains
66 
67 ! --------------------------------------------------------------------------------------------------
68 
69 subroutine geos_create(self, geom, conf)
70 
71 implicit none
72 type(geos_model), intent(inout) :: self
73 type(fv3jedi_geom), intent(in) :: geom
74 type(fckit_configuration), intent(in) :: conf
75 
76 type(esmf_vm) :: vm
77 type(mapl_capoptions) :: cap_options
78 type(mapl_communicators) :: mapl_comm
79 
80 integer :: rc, subcommunicator
81 
82 type(duration) :: dtstep
83 integer :: geos_dt, jedi_dt
84 character(len=20) :: ststep
85 character(len=:), allocatable :: str
86 logical :: esmf_logging
87 type(esmf_logkind_flag) :: esmf_logkind
88 
89 ! Grid dimensions
90 ! ---------------
91 self%isc = geom%isc
92 self%iec = geom%iec
93 self%jsc = geom%jsc
94 self%jec = geom%jec
95 self%npz = geom%npz
96 self%f_comm = geom%f_comm
97 
98 ! Set up options for the cap
99 ! --------------------------
100 cap_options = mapl_capoptions(cap_rc_file='CAP.rc')
101 cap_options%use_comm_world = .false.
102 cap_options%comm = geom%f_comm%communicator()
103 cap_options%npes_model = geom%f_comm%size()
104 
105 ! Create the MAPL_Cap object
106 ! --------------------------
107 self%Cap = mapl_cap( name='GCM', set_services=geos_gcm, cap_options = cap_options)
108 
109 ! MPI
110 ! ---
111 call self%cap%initialize_mpi()
112 
113 ! CapCom
114 ! ------
115 mapl_comm%mapl%comm = geom%f_comm%communicator()
116 mapl_comm%global%comm = geom%f_comm%communicator()
117 mapl_comm%esmf%comm = geom%f_comm%communicator()
118 mapl_comm%io%comm = mpi_comm_null
119 
120 call fill_comm(mapl_comm%mapl)
121 call fill_comm(mapl_comm%global)
122 call fill_comm(mapl_comm%esmf)
123 call fill_comm(mapl_comm%io)
124 
125 ! IO Server communicator (not on by default)
126 ! -----------------------------------------
127 subcommunicator = self%cap%create_member_subcommunicator(self%cap%get_comm_world(), rc=rc)
128 call self%cap%initialize_io_clients_servers(subcommunicator, rc = rc)
129 
130 if (conf%has('ESMF_Logging')) call conf%get_or_die('ESMF_Logging',esmf_logging)
131 ! Initialize ESMF
132 ! ---------------
133 if (esmf_logging) then
134  esmf_logkind = esmf_logkind_multi
135 else
136  esmf_logkind = esmf_logkind_none
137 end if
138 call esmf_initialize (vm=vm, logkindflag=esmf_logkind, mpicommunicator=mapl_comm%esmf%comm, rc=rc)
139 
140 ! Intialize the Cap GridComp
141 ! --------------------------
142 call self%cap%initialize_cap_gc(mapl_comm)
143 
144 ! GEOS SetServices
145 ! ----------------
146 call self%cap%cap_gc%set_services(rc = rc);
147 if (rc.ne.0) call abor1_ftn("geos_mod: set_services failed")
148 
149 ! GEOS Initialize
150 ! ---------------
151 call self%cap%cap_gc%initialize(rc = rc);
152 if (rc.ne.0) call abor1_ftn("geos_mod: initialize failed")
153 
154 ! Make sure required exports are allocated
155 ! ----------------------------------------
156 call allocate_exports(self)
157 
158 ! Time step checks and get number of GEOSsubsteps if any
159 ! ------------------------------------------------------
160 call conf%get_or_die("tstep",str)
161 ststep = str
162 deallocate(str)
163 dtstep = trim(ststep)
164 jedi_dt = int(duration_seconds(dtstep))
165 
166 geos_dt = self%cap%cap_gc%get_heartbeat_dt(rc = rc)
167 if (rc.ne.0) call abor1_ftn("geos_mod: error retrieving heartbeat")
168 
169 if (jedi_dt < geos_dt) then
170  call abor1_ftn("geos_mod: JEDI model time step should not be less than GEOS time step")
171 elseif (mod(jedi_dt,geos_dt) .ne. 0) then
172  call abor1_ftn("geos_mod: JEDI time step needs to be divisible by GEOS time step")
173 endif
174 
175 self%GEOSsubsteps = jedi_dt/geos_dt
176 
177 if (geom%f_comm%rank() == 0) then
178  print*, "There are ", int(self%GEOSsubsteps,2), " time steps of GEOS for each time step of JEDI"
179 endif
180 
181 self%start_time = self%cap%cap_gc%get_current_time(rc=rc)
182 self%saved_state = .false.
183 if (conf%has('reforecast')) call conf%get_or_die('reforecast',self%reforecast)
184 
185 end subroutine geos_create
186 
187 ! --------------------------------------------------------------------------------------------------
188 
189 subroutine geos_delete(self)
190 
191 implicit none
192 type(geos_model), intent(inout) :: self
193 
194 integer :: rc
195 
196 ! Finalize GEOS
197 ! -------------
198 call self%cap%cap_gc%destroy_state(rc = rc)
199 call self%cap%cap_gc%finalize(rc = rc);
200 if (rc.ne.0) call abor1_ftn("geos_mod: finalize failed")
201 call self%cap%finalize_io_clients_servers()
202 call mapl_finalize(comm=self%f_comm%communicator())
203 
204 ! Finalize ESMF
205 ! -------------
206 call esmf_finalize(endflag=esmf_end_keepmpi,rc=rc)
207 
208 end subroutine geos_delete
209 
210 ! --------------------------------------------------------------------------------------------------
211 
212 subroutine geos_initialize(self, state)
213 
214 implicit none
215 type(geos_model), target :: self
216 type(fv3jedi_state) :: state
217 integer :: rc
218 !There is no GEOS equivalent to JEDI initialize
219 
220 if (self%saved_state) then
221  call self%cap%rewind_model(self%start_time,rc=rc)
222  call self%cap%cap_gc%refresh_state(rc=rc)
223 end if
224 
225 if (self%reforecast .and. (.not.self%saved_state) ) then
226  call self%cap%cap_gc%record_state(rc=rc)
227  self%saved_state = .true.
228 end if
229 
230 
231 
232 end subroutine geos_initialize
233 
234 ! --------------------------------------------------------------------------------------------------
235 
236 subroutine geos_step(self, state, vdate)
237 
238 implicit none
239 type(geos_model), intent(inout) :: self
240 type(fv3jedi_state), intent(inout) :: state
241 type(datetime), intent(in) :: vdate !< Valid datetime after step
242 
243 integer :: n, rc
244 
245 !Convert JEDI state to GEOS state
246 !call state_to_geos( state, self )
247 
248 !Cycle GEOS through this time step of JEDI
249 do n = 1,self%GEOSsubsteps
250  rc = esmf_success
251  call self%cap%step_model(rc = rc)
252  if (rc.ne.0) call abor1_ftn("geos_mod: step_model failed")
253 enddo
254 
255 !Retrieve GEOS state and put into JEDI state
256 call geos_to_state( self, state )
257 
258 end subroutine geos_step
259 
260 ! --------------------------------------------------------------------------------------------------
261 
262 subroutine geos_finalize(self, state)
263 
264 implicit none
265 type(geos_model), target :: self
266 type(fv3jedi_state) :: state
267 
268 !There is no GEOS equivalent to JEDI finalize
269 
270 end subroutine geos_finalize
271 
272 ! --------------------------------------------------------------------------------------------------
273 
274 ! subroutine state_to_geos( state, self )
275 !
276 ! implicit none
277 ! type(fv3jedi_state), intent(in) :: state
278 ! type(geos_model), intent(inout) :: self
279 !
280 ! !TODO
281 !
282 ! ! 1. Overwrite the MKIAU exports?
283 ! ! 2. Dont compile with MKIAU and intercept EXT data?
284 ! ! 3. Cinderella grid comp for data?
285 !
286 ! end subroutine state_to_geos
287 
288 ! --------------------------------------------------------------------------------------------------
289 
290 subroutine geos_to_state( self, state )
291 
292 implicit none
293 type(geos_model), intent(in) :: self
294 type(fv3jedi_state), intent(inout) :: state
295 
296 integer :: num_items, i, rc, rank, lb(3), ub(3), fnpz
297 type(esmf_field) :: field
298 character(len=ESMF_MAXSTR), allocatable :: item_names(:)
299 real(kind=esmf_kind_r4), pointer :: farrayptr2(:,:)
300 real(kind=esmf_kind_r4), pointer :: farrayptr3(:,:,:)
301 character(len=field_clen) :: fv3jedi_name
302 type(fv3jedi_field), pointer :: field_ptr
303 real(kind=kind_real), allocatable, dimension(:,:,:) :: field_geos
304 
305 
306 ! Array to hold output from GEOS in JEDI precision
307 ! ------------------------------------------------
308 allocate(field_geos(self%isc:self%iec, self%jsc:self%jec, self%npz+1))
309 
310 
311 ! Get number of items
312 ! -------------------
313 call esmf_stateget(self%cap%cap_gc%export_state, itemcount = num_items, rc = rc)
314 if (rc.ne.0) call abor1_ftn("geos_to_state: ESMF_StateGet itemcount failed")
315 
316 
317 ! Get names of the items
318 ! ----------------------
319 allocate(item_names(num_items))
320 call esmf_stateget(self%cap%cap_gc%export_state, itemnamelist = item_names, rc = rc)
321 if (rc.ne.0) call abor1_ftn("geos_to_state: ESMF_StateGet itemnamelist failed")
322 
323 
324 ! Loop over states coming from GEOS and convert to JEDI state
325 ! -----------------------------------------------------------
326 do i = 1, num_items
327 
328  ! Create map between GEOS name and fv3-jedi name
329  ! ----------------------------------------------
330  select case (trim(item_names(i)))
331 
332  ! The below need to be listed in Cap.rc as e.g.
333  ! CAP_EXPORTS:
334  ! U_DGRID,DYN
335  ! Q,MOIST
336  ! ::
337 
338  ! DYN
339  case ("U_DGRID")
340  fv3jedi_name = 'ud'
341  case ("V_DGRID")
342  fv3jedi_name = 'vd'
343  case ("PT")
344  fv3jedi_name = 'pt'
345  case ("PE")
346  fv3jedi_name = 'pe'
347  case ("T")
348  fv3jedi_name = 't'
349  case ("DELP")
350  fv3jedi_name = 'delp'
351  case ("PS")
352  fv3jedi_name = 'ps'
353  case ("TA")
354  fv3jedi_name = 'ts'
355  case ("DZ")
356  fv3jedi_name = 'delz'
357  case ("W")
358  fv3jedi_name = 'w'
359 
360  ! AGCM
361  case ("PHIS")
362  fv3jedi_name = 'phis'
363  case ("QITOT")
364  fv3jedi_name = 'ice_wat'
365  case ("QLTOT")
366  fv3jedi_name = 'liq_wat'
367  case ("VARFLT")
368  fv3jedi_name = 'varflt'
369 
370  ! MOIST
371  case ("Q")
372  fv3jedi_name = 'sphum'
373  case ("QICN")
374  fv3jedi_name = "qicn"
375  case ("QLCN")
376  fv3jedi_name = "qlcn"
377  case ("QILS")
378  fv3jedi_name = "qils"
379  case ("QLLS")
380  fv3jedi_name = "qlls"
381  case ("QSTOT")
382  fv3jedi_name = "qs"
383  case ("QRTOT")
384  fv3jedi_name = "qr"
385  case ("QCLSX0")
386  fv3jedi_name = 'qls'
387  case ("QCCNX0")
388  fv3jedi_name = 'qcn'
389  case ("CLCNX0")
390  fv3jedi_name = 'cfcn'
391  case ("KCBL_moist")
392  fv3jedi_name = 'kcbl'
393  case ("TS_moist")
394  fv3jedi_name = 'tsm'
395  case ("KHl_moist")
396  fv3jedi_name = 'khl'
397  case ("KHu_moist")
398  fv3jedi_name = 'khu'
399 
400  ! CHEMISTRY
401  case ("O3")
402  fv3jedi_name = 'o3ppmv'
403 
404  ! SURFACE
405  case ("FRLAND")
406  fv3jedi_name = 'frland'
407  case ("FRLANDICE")
408  fv3jedi_name = 'frlandice'
409  case ("FRLAKE")
410  fv3jedi_name = 'frlake'
411  case ("FROCEAN")
412  fv3jedi_name = 'frocean'
413  case ("FRACI")
414  fv3jedi_name = 'frseaice'
415  case ("USTAR")
416  fv3jedi_name = 'ustar'
417  case ("BSTAR")
418  fv3jedi_name = 'bstar'
419  case ("CM")
420  fv3jedi_name = 'cm'
421  case ("CT")
422  fv3jedi_name = 'ct'
423  case ("CQ")
424  fv3jedi_name = 'cq'
425  case ("U10N")
426  fv3jedi_name = 'u_srf'
427  case ("V10N")
428  fv3jedi_name = 'v_srf'
429  case ('SNOMAS')
430  fv3jedi_name = 'sheleg'
431  case ('TSOIL1')
432  fv3jedi_name = 'soilt'
433  case ('WET1')
434  fv3jedi_name = 'soilm'
435 
436  !TURBULENCE
437  case ('ZPBL')
438  fv3jedi_name = 'zpbl'
439 
440  ! NO MAP
441  case default
442  call abor1_ftn("geos_to_state no map between GEOS name and JEDI name "//trim(item_names(i)))
443 
444  end select
445 
446 
447  ! Only need to extract field from GEOS if fv3-jedi needs it
448  ! ---------------------------------------------------------
449  if (state%has_field(trim(fv3jedi_name))) then
450 
451  !Get field from the state
452  call esmf_stateget(self%cap%cap_gc%export_state, item_names(i), field, rc = rc)
453  if (rc.ne.0) call abor1_ftn("geos_to_state: ESMF_StateGet field failed")
454 
455  !Validate the field
456  call esmf_fieldvalidate(field, rc = rc)
457  if (rc.ne.0) call abor1_ftn("geos_to_state: ESMF_FieldValidate failed")
458 
459  !Get the field rank
460  call esmf_fieldget(field, rank = rank, rc = rc)
461  if (rc.ne.0) call abor1_ftn("geos_to_state: ESMF_FieldGet rank failed")
462 
463  !Convert field to pointer and pointer bounds
464  field_geos = 0.0_kind_real
465  if (rank == 2) then
466 
467  call esmf_fieldget( field, 0, farrayptr = farrayptr2, totallbound = lb(1:2), totalubound = ub(1:2), rc = rc )
468  if (rc.ne.0) call abor1_ftn("geos_to_state: ESMF_FieldGet 2D failed")
469 
470  fnpz = 1
471  field_geos(self%isc:self%iec,self%jsc:self%jec,1) = farrayptr2(lb(1):ub(1),lb(2):ub(2))
472  nullify(farrayptr2)
473 
474  elseif (rank == 3) then
475 
476  call esmf_fieldget( field, 0, farrayptr = farrayptr3, totallbound = lb, totalubound = ub, rc = rc )
477  if (rc.ne.0) call abor1_ftn("geos_to_state: ESMF_FieldGet 3D failed")
478 
479  fnpz = ub(3)-lb(3)+1
480  field_geos(self%isc:self%iec,self%jsc:self%jec,1:fnpz) = farrayptr3(lb(1):ub(1),lb(2):ub(2),lb(3):ub(3))
481  nullify(farrayptr3)
482 
483  else
484 
485  call abor1_ftn("geos_mod: can only handle rank 2 or rank 3 fields from GEOS")
486 
487  endif
488 
489  ! Check that dimensions match
490  if ((ub(1)-lb(1)+1 .ne. self%iec-self%isc+1) .or. (ub(2)-lb(2)+1 .ne. self%jec-self%jsc+1) ) then
491  call abor1_ftn("geos_to_state: dimension mismatch between JEDI and GEOS horizontal grid")
492  endif
493 
494  ! Get pointer to fv3-jedi side field
495  call state%get_field(trim(fv3jedi_name), field_ptr)
496 
497  if (field_ptr%npz .ne. fnpz) &
498  call abor1_ftn("geos_to_state: dimension mismatch between JEDI and GEOS vertical grid")
499 
500  ! Copy from GEOS to fv3-jedi
501  field_ptr%array(:,:,1:fnpz) = field_geos(:,:,1:fnpz)
502  endif
503 
504 end do
505 
506 deallocate(item_names)
507 
508 end subroutine geos_to_state
509 
510 ! --------------------------------------------------------------------------------------------------
511 
512 subroutine allocate_exports(self)
513 
514 implicit none
515 type(geos_model), intent(inout) :: self
516 
517 integer :: num_items, i, rc
518 type(esmf_field) :: field
519 character(len=ESMF_MAXSTR), allocatable :: item_names(:)
520 
521 
522  ! MAPL exports are only allocated if a connection is found between the
523  ! export and something else, e.g. HISTORY. This routine creates this
524  ! link between exports and MAPL_CapGridComp
525 
526  ! Get number of items
527  ! -------------------
528  call esmf_stateget(self%cap%cap_gc%export_state, itemcount = num_items, rc = rc)
529  if (rc.ne.0) call abor1_ftn("allocate_exports: ESMF_StateGet itemcount failed")
530 
531 
532  ! Get names of the items
533  ! ----------------------
534  allocate(item_names(num_items))
535  call esmf_stateget(self%cap%cap_gc%export_state, itemnamelist = item_names, rc = rc)
536  if (rc.ne.0) call abor1_ftn("allocate_exports: ESMF_StateGet items failed")
537 
538 
539  ! Loop over states and allocate coupling
540  ! --------------------------------------
541  do i = 1, num_items
542 
543  !Get field from the state
544  call esmf_stateget(self%cap%cap_gc%export_state, item_names(i), field, rc = rc)
545  if (rc.ne.0) call abor1_ftn("allocate_exports: ESMF_StateGet field failed")
546 
547  !Validate the field
548  call esmf_fieldvalidate(field, rc = rc)
549  if (rc.ne.0) call abor1_ftn("allocate_exports: ESMF_FieldValidate failed")
550 
551  call mapl_allocatecoupling(field, rc)
552  if (rc.ne.0) call abor1_ftn("allocate_exports: MAPL_AllocateCoupling failed")
553 
554  end do
555 
556  deallocate(item_names)
557 
558 end subroutine allocate_exports
559 
560 ! --------------------------------------------------------------------------------------------------
561 
562 subroutine fill_comm(mcomm)
563 
564  type (MAPL_Communicator), intent(inout) :: mcomm
565 
566  integer :: comm
567  integer :: status
568 
569  comm = mcomm%comm
570  if (comm == mpi_comm_null) then
571  mcomm%size = 0
572  mcomm%rank = mpi_undefined
573  mcomm%root = mpi_undefined
574  else
575  call mpi_comm_size(comm, mcomm%size,status)
576  _verify(status)
577  call mpi_comm_rank(comm, mcomm%rank,status)
578  _verify(status)
579  mcomm%root = 0
580  end if
581 end subroutine fill_comm
582 
583 ! --------------------------------------------------------------------------------------------------
584 
585 
586 end module fv3jedi_geos_mod
fv3jedi_geos_mod
Definition: fv3jedi_geos_mod.F90:10
fv3jedi_state_mod::fv3jedi_state
Fortran derived type to hold FV3JEDI state.
Definition: fv3jedi_state_mod.F90:30
fv3jedi_field_mod
Definition: fv3jedi_field_mod.f90:6
fv3jedi_state_mod
Definition: fv3jedi_state_mod.F90:6
fv3jedi_geom_mod
Fortran module handling geometry for the FV3 model.
Definition: fv3jedi_geom_mod.f90:8
fv3jedi_geos_mod::geos_initialize
subroutine, public geos_initialize(self, state)
Definition: fv3jedi_geos_mod.F90:213
fv3jedi_geos_mod::fill_comm
subroutine fill_comm(mcomm)
Definition: fv3jedi_geos_mod.F90:563
fv3jedi_increment_mod
Definition: fv3jedi_increment_mod.F90:6
fv3jedi_geom_mod::fv3jedi_geom
Fortran derived type to hold geometry data for the FV3JEDI model.
Definition: fv3jedi_geom_mod.f90:46
fv3jedi_geos_mod::geos_create
subroutine, public geos_create(self, geom, conf)
Definition: fv3jedi_geos_mod.F90:70
fv3jedi_geos_mod::field_attributes
Definition: fv3jedi_geos_mod.F90:60
fv3jedi_kinds_mod::kind_real
integer, parameter, public kind_real
Definition: fv3jedi_kinds_mod.f90:14
fv3jedi_geos_mod::geos_finalize
subroutine, public geos_finalize(self, state)
Definition: fv3jedi_geos_mod.F90:263
fv3jedi_geos_mod::geos_to_state
subroutine geos_to_state(self, state)
Definition: fv3jedi_geos_mod.F90:291
fv3jedi_field_mod::fv3jedi_field
Definition: fv3jedi_field_mod.f90:36
fv3jedi_kinds_mod
Definition: fv3jedi_kinds_mod.f90:6
fv3jedi_geos_mod::geos_step
subroutine, public geos_step(self, state, vdate)
Definition: fv3jedi_geos_mod.F90:237
fv3jedi_geos_mod::allocate_exports
subroutine allocate_exports(self)
Definition: fv3jedi_geos_mod.F90:513
fv3jedi_geos_mod::geos_model
Fortran derived type to hold model definition.
Definition: fv3jedi_geos_mod.F90:48
fv3jedi_geos_mod::geos_delete
subroutine, public geos_delete(self)
Definition: fv3jedi_geos_mod.F90:190
fv3jedi_increment_mod::fv3jedi_increment
Definition: fv3jedi_increment_mod.F90:34
fv3jedi_field_mod::field_clen
integer, parameter, public field_clen
Definition: fv3jedi_field_mod.f90:31