From f59b5c3d6846787f5c72938a83b77b78fa032dea Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Thu, 28 Jan 2021 11:41:28 -0700 Subject: [PATCH 1/4] changes for creating masks at runtime --- cicecore/cicedynB/general/ice_init.F90 | 2 +- cicecore/cicedynB/infrastructure/ice_grid.F90 | 285 ------ cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 | 844 +++++++++--------- .../drivers/nuopc/cmeps/ice_comp_nuopc.F90 | 315 ++----- 4 files changed, 484 insertions(+), 962 deletions(-) diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index bde2829da..3a5a8ee03 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -1604,7 +1604,7 @@ subroutine input_data grid_type /= 'rectangular' .and. & grid_type /= 'cpom_grid' .and. & grid_type /= 'regional' .and. & - grid_type /= 'latlon' ) then + grid_type /= 'setmask' ) then if (my_task == master_task) write(nu_diag,*) subname//' ERROR: unknown grid_type=',trim(grid_type) abort_list = trim(abort_list)//":20" endif diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index 34b37cf29..ce1af48d8 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -368,11 +368,6 @@ subroutine init_grid2 else call popgrid ! read POP grid lengths directly endif -#ifdef CESMCOUPLED - elseif (trim(grid_type) == 'latlon') then - call latlongrid ! lat lon grid for sequential CESM (CAM mode) - return -#endif elseif (trim(grid_type) == 'cpom_grid') then call cpomgrid ! cpom model orca1 type grid else @@ -882,286 +877,6 @@ subroutine popgrid_nc end subroutine popgrid_nc -#ifdef CESMCOUPLED -!======================================================================= - -! Read in kmt file that matches CAM lat-lon grid and has single column -! functionality -! author: Mariana Vertenstein -! 2007: Elizabeth Hunke upgraded to netcdf90 and cice ncdf calls - - subroutine latlongrid - -! use ice_boundary - use ice_domain_size - use ice_scam, only : scmlat, scmlon, single_column - use ice_constants, only: c0, c1, p5, p25, & - field_loc_center, field_type_scalar, radius -#ifdef USE_NETCDF - use netcdf -#endif - - integer (kind=int_kind) :: & - i, j, iblk - - integer (kind=int_kind) :: & - ni, nj, ncid, dimid, varid, ier - - type (block) :: & - this_block ! block information for current block - - integer (kind=int_kind) :: & - ilo,ihi,jlo,jhi ! beginning and end of physical domain - - real (kind=dbl_kind) :: & - closelat, & ! Single-column latitude value - closelon, & ! Single-column longitude value - closelatidx, & ! Single-column latitude index to retrieve - closelonidx ! Single-column longitude index to retrieve - - integer (kind=int_kind) :: & - start(2), & ! Start index to read in - count(2) ! Number of points to read in - - integer (kind=int_kind) :: & - start3(3), & ! Start index to read in - count3(3) ! Number of points to read in - - integer (kind=int_kind) :: & - status ! status flag - - real (kind=dbl_kind), allocatable :: & - lats(:),lons(:),pos_lons(:), glob_grid(:,:) ! temporaries - - real (kind=dbl_kind) :: & - pos_scmlon,& ! temporary - pi, & - puny, & - scamdata ! temporary - - character(len=*), parameter :: subname = '(lonlatgrid)' - -#ifdef USE_NETCDF - !----------------------------------------------------------------- - ! - kmt file is actually clm fractional land file - ! - Determine consistency of dimensions - ! - Read in lon/lat centers in degrees from kmt file - ! - Read in ocean from "kmt" file (1 for ocean, 0 for land) - !----------------------------------------------------------------- - - call icepack_query_parameters(pi_out=pi, puny_out=puny) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & - file=__FILE__, line=__LINE__) - - ! Determine dimension of domain file and check for consistency - - if (my_task == master_task) then - call ice_open_nc(kmt_file, ncid) - - status = nf90_inq_dimid (ncid, 'ni', dimid) - status = nf90_inquire_dimension(ncid, dimid, len=ni) - status = nf90_inq_dimid (ncid, 'nj', dimid) - status = nf90_inquire_dimension(ncid, dimid, len=nj) - end if - - ! Determine start/count to read in for either single column or global lat-lon grid - ! If single_column, then assume that only master_task is used since there is only one task - - if (single_column) then - ! Check for consistency - if (my_task == master_task) then - if ((nx_global /= 1).or. (ny_global /= 1)) then - write(nu_diag,*) 'Because you have selected the column model flag' - write(nu_diag,*) 'Please set nx_global=ny_global=1 in file' - write(nu_diag,*) 'ice_domain_size.F and recompile' - call abort_ice (subname//'ERROR: check nx_global, ny_global') - endif - end if - - ! Read in domain file for single column - allocate(lats(nj)) - allocate(lons(ni)) - allocate(pos_lons(ni)) - allocate(glob_grid(ni,nj)) - - start3=(/1,1,1/) - count3=(/ni,nj,1/) - status = nf90_inq_varid(ncid, 'xc' , varid) - if (status /= nf90_noerr) call abort_ice (subname//' inq_varid xc') - status = nf90_get_var(ncid, varid, glob_grid, start3, count3) - if (status /= nf90_noerr) call abort_ice (subname//' get_var xc') - do i = 1,ni - lons(i) = glob_grid(i,1) - end do - - status = nf90_inq_varid(ncid, 'yc' , varid) - if (status /= nf90_noerr) call abort_ice (subname//' inq_varid yc') - status = nf90_get_var(ncid, varid, glob_grid, start3, count3) - if (status /= nf90_noerr) call abort_ice (subname//' get_var yc') - do j = 1,nj - lats(j) = glob_grid(1,j) - end do - - ! convert lons array and scmlon to 0,360 and find index of value closest to 0 - ! and obtain single-column longitude/latitude indices to retrieve - - pos_lons(:)= mod(lons(:) + 360._dbl_kind,360._dbl_kind) - pos_scmlon = mod(scmlon + 360._dbl_kind,360._dbl_kind) - start(1) = (MINLOC(abs(pos_lons-pos_scmlon),dim=1)) - start(2) = (MINLOC(abs(lats -scmlat ),dim=1)) - - deallocate(lats) - deallocate(lons) - deallocate(pos_lons) - deallocate(glob_grid) - - status = nf90_inq_varid(ncid, 'xc' , varid) - if (status /= nf90_noerr) call abort_ice (subname//' inq_varid xc') - status = nf90_get_var(ncid, varid, scamdata, start) - if (status /= nf90_noerr) call abort_ice (subname//' get_var xc') - TLON = scamdata - status = nf90_inq_varid(ncid, 'yc' , varid) - if (status /= nf90_noerr) call abort_ice (subname//' inq_varid yc') - status = nf90_get_var(ncid, varid, scamdata, start) - if (status /= nf90_noerr) call abort_ice (subname//' get_var yc') - TLAT = scamdata - status = nf90_inq_varid(ncid, 'area' , varid) - if (status /= nf90_noerr) call abort_ice (subname//' inq_varid area') - status = nf90_get_var(ncid, varid, scamdata, start) - if (status /= nf90_noerr) call abort_ice (subname//' get_var are') - tarea = scamdata - status = nf90_inq_varid(ncid, 'mask' , varid) - if (status /= nf90_noerr) call abort_ice (subname//' inq_varid mask') - status = nf90_get_var(ncid, varid, scamdata, start) - if (status /= nf90_noerr) call abort_ice (subname//' get_var mask') - hm = scamdata - status = nf90_inq_varid(ncid, 'frac' , varid) - if (status /= nf90_noerr) call abort_ice (subname//' inq_varid frac') - status = nf90_get_var(ncid, varid, scamdata, start) - if (status /= nf90_noerr) call abort_ice (subname//' get_var frac') - ocn_gridcell_frac = scamdata - else - ! Check for consistency - if (my_task == master_task) then - if (nx_global /= ni .and. ny_global /= nj) then - write(nu_diag,*) 'latlongrid: ni,nj = ',ni,nj - write(nu_diag,*) 'latlongrid: nx_g,ny_g = ',nx_global, ny_global - call abort_ice (subname//'ERROR: ni,nj not equal to nx_global,ny_global') - end if - end if - - ! Read in domain file for global lat-lon grid - call ice_read_nc(ncid, 1, 'xc' , TLON , diag=.true.) - call ice_read_nc(ncid, 1, 'yc' , TLAT , diag=.true.) - call ice_read_nc(ncid, 1, 'area', tarea , diag=.true., & - field_loc=field_loc_center,field_type=field_type_scalar) - call ice_read_nc(ncid, 1, 'mask', hm , diag=.true.) - call ice_read_nc(ncid, 1, 'frac', ocn_gridcell_frac, diag=.true.) - end if - - if (my_task == master_task) then - call ice_close_nc(ncid) - end if - - !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,i,j) - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - - do j = jlo, jhi - do i = ilo, ihi - ! Convert from degrees to radians - TLON(i,j,iblk) = pi*TLON(i,j,iblk)/180._dbl_kind - - ! Convert from degrees to radians - TLAT(i,j,iblk) = pi*TLAT(i,j,iblk)/180._dbl_kind - - ! Convert from radians^2 to m^2 - ! (area in domain file is in radians^2 and tarea is in m^2) - tarea(i,j,iblk) = tarea(i,j,iblk) * (radius*radius) - end do - end do - end do - !$OMP END PARALLEL DO - - !----------------------------------------------------------------- - ! Calculate various geometric 2d arrays - ! The U grid (velocity) is not used when run with sequential CAM - ! because we only use thermodynamic sea ice. However, ULAT is used - ! in the default initialization of CICE so we calculate it here as - ! a "dummy" so that CICE will initialize with ice. If a no ice - ! initialization is OK (or desired) this can be commented out and - ! ULAT will remain 0 as specified above. ULAT is located at the - ! NE corner of the grid cell, TLAT at the center, so here ULAT is - ! hacked by adding half the latitudinal spacing (in radians) to - ! TLAT. - !----------------------------------------------------------------- - - !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,i,j) - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - - do j = jlo, jhi - do i = ilo, ihi - - if (ny_global == 1) then - uarea(i,j,iblk) = tarea(i,j, iblk) - else - uarea(i,j,iblk) = p25* & - (tarea(i,j, iblk) + tarea(i+1,j, iblk) & - + tarea(i,j+1,iblk) + tarea(i+1,j+1,iblk)) - endif - tarear(i,j,iblk) = c1/tarea(i,j,iblk) - uarear(i,j,iblk) = c1/uarea(i,j,iblk) - tinyarea(i,j,iblk) = puny*tarea(i,j,iblk) - - if (single_column) then - ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/nj) - else - if (ny_global == 1) then - ULAT (i,j,iblk) = TLAT(i,j,iblk) - else - ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/ny_global) - endif - endif - ULON (i,j,iblk) = c0 - ANGLE (i,j,iblk) = c0 - - ANGLET(i,j,iblk) = c0 - HTN (i,j,iblk) = 1.e36_dbl_kind - HTE (i,j,iblk) = 1.e36_dbl_kind - dxt (i,j,iblk) = 1.e36_dbl_kind - dyt (i,j,iblk) = 1.e36_dbl_kind - dxu (i,j,iblk) = 1.e36_dbl_kind - dyu (i,j,iblk) = 1.e36_dbl_kind - dxhy (i,j,iblk) = 1.e36_dbl_kind - dyhx (i,j,iblk) = 1.e36_dbl_kind - cyp (i,j,iblk) = 1.e36_dbl_kind - cxp (i,j,iblk) = 1.e36_dbl_kind - cym (i,j,iblk) = 1.e36_dbl_kind - cxm (i,j,iblk) = 1.e36_dbl_kind - enddo - enddo - enddo - !$OMP END PARALLEL DO - - call makemask -#else - call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', & - file=__FILE__, line=__LINE__) -#endif - - end subroutine latlongrid -#endif - !======================================================================= ! Regular rectangular grid and mask diff --git a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 index 917774908..a585dd331 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 @@ -1,444 +1,446 @@ -!======================================================================= -! -! This module contains the CICE initialization routine that sets model -! parameters and initializes the grid and CICE state variables. -! -! authors Elizabeth C. Hunke, LANL -! William H. Lipscomb, LANL -! Philip W. Jones, LANL -! -! 2006: Converted to free form source (F90) by Elizabeth Hunke -! 2008: E. Hunke moved ESMF code to its own driver - - module CICE_InitMod - - use ice_kinds_mod - use ice_exit, only: abort_ice - use ice_fileunits, only: init_fileunits, nu_diag - use icepack_intfc, only: icepack_aggregate - use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist - use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave - use icepack_intfc, only: icepack_configure - use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted - use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags, & - icepack_query_tracer_indices, icepack_query_tracer_sizes - - implicit none - private - public :: cice_init +module CICE_InitMod -!======================================================================= + ! Initialize CICE model. - contains + use ice_kinds_mod + use ice_exit, only: abort_ice + use ice_fileunits, only: init_fileunits, nu_diag + use icepack_intfc, only: icepack_aggregate + use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist + use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave + use icepack_intfc, only: icepack_configure + use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted + use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags + use icepack_intfc, only: icepack_query_tracer_indices, icepack_query_tracer_sizes + implicit none + private + public :: cice_init + +!======================================================================= +contains !======================================================================= -! -! Initialize CICE model. - - subroutine cice_init - - ! Initialize the basic state, grid and all necessary parameters for - ! running the CICE model. - - use ice_arrays_column, only: hin_max, c_hi_range, alloc_arrays_column - use ice_arrays_column, only: floe_rad_l, floe_rad_c, & - floe_binwidth, c_fsd_range - use ice_state, only: alloc_state - use ice_flux_bgc, only: alloc_flux_bgc - use ice_calendar, only: dt, dt_dyn, time, istep, istep1, write_ic, & - init_calendar, calendar - use ice_communicate, only: my_task, master_task - use ice_diagnostics, only: init_diags - use ice_domain, only: init_domain_blocks - use ice_domain_size, only: ncat, nfsd - use ice_dyn_eap, only: init_eap, alloc_dyn_eap - use ice_dyn_shared, only: kdyn, init_evp, alloc_dyn_shared - use ice_flux, only: init_coupler_flux, init_history_therm, & - init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux - use ice_forcing, only: init_forcing_ocn - use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, & - faero_default, faero_optics, alloc_forcing_bgc, fiso_default - use ice_grid, only: init_grid1, init_grid2, alloc_grid - use ice_history, only: init_hist, accum_hist - use ice_restart_shared, only: restart, runtype - use ice_init, only: input_data, init_state - use ice_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc, input_zbgc, count_tracers - use ice_kinds_mod - use ice_restoring, only: ice_HaloRestore_init - use ice_timers, only: timer_total, init_ice_timers, ice_timer_start - use ice_transport_driver, only: init_transport - - logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers, & - tr_iso, tr_fsd, wave_spec - character(len=*), parameter :: subname = '(cice_init)' - - call init_fileunits ! unit numbers - - call icepack_configure() ! initialize icepack - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(trim(subname), & - file=__FILE__,line= __LINE__) - - call input_data ! namelist variables - call input_zbgc ! vertical biogeochemistry namelist - call count_tracers ! count tracers - - call init_domain_blocks ! set up block decomposition - call init_grid1 ! domain distribution - call alloc_grid ! allocate grid arrays - call alloc_arrays_column ! allocate column arrays - call alloc_state ! allocate state arrays - call alloc_dyn_shared ! allocate dyn shared arrays - call alloc_flux_bgc ! allocate flux_bgc arrays - call alloc_flux ! allocate flux arrays - call init_ice_timers ! initialize all timers - call ice_timer_start(timer_total) ! start timing entire run - call init_grid2 ! grid variables - call init_zbgc ! vertical biogeochemistry initialization - call init_calendar ! initialize some calendar stuff - call init_hist (dt) ! initialize output history file - - if (kdyn == 2) then - call alloc_dyn_eap ! allocate dyn_eap arrays - call init_eap (dt_dyn) ! define eap dynamics parameters, variables - else ! for both kdyn = 0 or 1 - call init_evp (dt_dyn) ! define evp dynamics parameters, variables - endif - - call init_coupler_flux ! initialize fluxes exchanged with coupler - call init_thermo_vertical ! initialize vertical thermodynamics - - call icepack_init_itd(ncat=ncat, hin_max=hin_max) ! ice thickness distribution - if (my_task == master_task) then - call icepack_init_itd_hist(ncat=ncat, hin_max=hin_max, c_hi_range=c_hi_range) ! output - endif - - call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(trim(subname), & - file=__FILE__,line= __LINE__) - - if (tr_fsd) call icepack_init_fsd_bounds (nfsd, & ! floe size distribution + + subroutine cice_init1() + + ! Initialize the basic state, grid and all necessary parameters for + ! running the CICE model. + + use ice_init , only: input_data + use ice_init_column , only: input_zbgc, count_tracers + use ice_grid , only: init_grid1, alloc_grid + use ice_domain , only: init_domain_blocks + use ice_arrays_column , only: alloc_arrays_column + use ice_state , only: alloc_state + use ice_dyn_shared , only: alloc_dyn_shared + use ice_flux_bgc , only: alloc_flux_bgc + use ice_flux , only: alloc_flux + use ice_timers , only: timer_total, init_ice_timers, ice_timer_start + + character(len=*), parameter :: subname = '(cice_init1)' + !---------------------------------------------------- + + call init_fileunits ! unit numbers + call icepack_configure() ! initialize icepack + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + call input_data ! namelist variables + call input_zbgc ! vertical biogeochemistry namelist + call count_tracers ! count tracers + + call init_domain_blocks ! set up block decomposition + call init_grid1 ! domain distribution + call alloc_grid ! allocate grid arrays + call alloc_arrays_column ! allocate column arrays + call alloc_state ! allocate state arrays + call alloc_dyn_shared ! allocate dyn shared arrays + call alloc_flux_bgc ! allocate flux_bgc arrays + call alloc_flux ! allocate flux arrays + call init_ice_timers ! initialize all timers + call ice_timer_start(timer_total) ! start timing entire run + + end subroutine cice_init1 + + !======================================================================= + subroutine cice_init2() + + ! Initialize the basic state, and all necessary parameters for + ! running the CICE model. + + use ice_arrays_column , only: hin_max, c_hi_range + use ice_arrays_column , only: floe_rad_l, floe_rad_c, floe_binwidth, c_fsd_range + use ice_calendar , only: dt, dt_dyn, time, istep, istep1, write_ic, init_calendar, calendar + use ice_communicate , only: my_task, master_task + use ice_diagnostics , only: init_diags + use ice_domain_size , only: ncat, nfsd + use ice_dyn_eap , only: init_eap, alloc_dyn_eap + use ice_dyn_shared , only: kdyn, init_evp + use ice_flux , only: init_coupler_flux, init_history_therm + use ice_flux , only: init_history_dyn, init_flux_atm, init_flux_ocn + use ice_forcing , only: init_forcing_ocn + use ice_forcing_bgc , only: get_forcing_bgc, get_atm_bgc + use ice_forcing_bgc , only: faero_default, faero_optics, alloc_forcing_bgc, fiso_default + use ice_history , only: init_hist, accum_hist + use ice_restart_shared , only: restart, runtype + use ice_init , only: input_data, init_state + use ice_init_column , only: init_thermo_vertical, init_shortwave, init_zbgc + use ice_restoring , only: ice_HaloRestore_init + use ice_timers , only: timer_total, init_ice_timers, ice_timer_start + use ice_transport_driver , only: init_transport + + logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers + logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec + character(len=*), parameter :: subname = '(cice_init2)' + !---------------------------------------------------- + + call init_zbgc ! vertical biogeochemistry initialization + call init_calendar ! initialize some calendar stuff + call init_hist (dt) ! initialize output history file + + if (kdyn == 2) then + call alloc_dyn_eap ! allocate dyn_eap arrays + call init_eap (dt_dyn) ! define eap dynamics parameters, variables + else ! for both kdyn = 0 or 1 + call init_evp (dt_dyn) ! define evp dynamics parameters, variables + endif + + call init_coupler_flux ! initialize fluxes exchanged with coupler + call init_thermo_vertical ! initialize vertical thermodynamics + + call icepack_init_itd(ncat=ncat, hin_max=hin_max) ! ice thickness distribution + if (my_task == master_task) then + call icepack_init_itd_hist(ncat=ncat, hin_max=hin_max, c_hi_range=c_hi_range) ! output + endif + + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (tr_fsd) call icepack_init_fsd_bounds (nfsd, & ! floe size distribution floe_rad_l, & ! fsd size lower bound in m (radius) floe_rad_c, & ! fsd size bin centre in m (radius) floe_binwidth, & ! fsd size bin width in m (radius) c_fsd_range, & ! string for history output write_diags=(my_task == master_task)) ! write diag on master only - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - call calendar(time) ! determine the initial date - - ! TODO: - why is this being called when you are using CMEPS? - call init_forcing_ocn(dt) ! initialize sss and sst from data - - call init_state ! initialize the ice state - call init_transport ! initialize horizontal transport - call ice_HaloRestore_init ! restored boundary conditions - - call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & - wave_spec_out=wave_spec) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(trim(subname), & - file=__FILE__,line= __LINE__) - - if (skl_bgc .or. z_tracers) call alloc_forcing_bgc ! allocate biogeochemistry arrays - - call init_restart ! initialize restart variables - call init_diags ! initialize diagnostic output points - call init_history_therm ! initialize thermo history variables - call init_history_dyn ! initialize dynamic history variables - - call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) - call icepack_query_tracer_flags(tr_iso_out=tr_iso) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(trim(subname), & - file=__FILE__,line= __LINE__) - - if (tr_aero .or. tr_zaero) then - call faero_optics !initialize aerosol optical property tables - end if - - ! Initialize shortwave components using swdn from previous timestep - ! if restarting. These components will be scaled to current forcing - ! in prep_radiation. - - if (trim(runtype) == 'continue' .or. restart) then - call init_shortwave ! initialize radiative transfer - end if - - !-------------------------------------------------------------------- - ! coupler communication or forcing data initialization - !-------------------------------------------------------------------- - - if (z_tracers) call get_atm_bgc ! biogeochemistry - - if (runtype == 'initial' .and. .not. restart) then - call init_shortwave ! initialize radiative transfer using current swdn - end if - - call init_flux_atm ! initialize atmosphere fluxes sent to coupler - call init_flux_ocn ! initialize ocean fluxes sent to coupler - - end subroutine cice_init - -!======================================================================= - - subroutine init_restart - - use ice_arrays_column, only: dhsn - use ice_blocks, only: nx_block, ny_block - use ice_calendar, only: time, calendar - use ice_constants, only: c0 - use ice_domain, only: nblocks - use ice_domain_size, only: ncat, n_iso, n_aero, nfsd - use ice_dyn_eap, only: read_restart_eap - use ice_dyn_shared, only: kdyn - use ice_grid, only: tmask - use ice_init, only: ice_ic - use ice_init_column, only: init_age, init_FY, init_lvl, & - init_meltponds_cesm, init_meltponds_lvl, init_meltponds_topo, & - init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd - use ice_restart_column, only: restart_age, read_restart_age, & - restart_FY, read_restart_FY, restart_lvl, read_restart_lvl, & - restart_pond_cesm, read_restart_pond_cesm, & - restart_pond_lvl, read_restart_pond_lvl, & - restart_pond_topo, read_restart_pond_topo, & - restart_fsd, read_restart_fsd, & - restart_iso, read_restart_iso, & - restart_aero, read_restart_aero, & - restart_hbrine, read_restart_hbrine, & - restart_zsal, restart_bgc - use ice_restart_driver, only: restartfile - use ice_restart_shared, only: runtype, restart - use ice_state ! almost everything - - integer(kind=int_kind) :: & + call calendar(time) ! determine the initial date + + ! TODO: - why is this being called when you are using CMEPS? + call init_forcing_ocn(dt) ! initialize sss and sst from data + + call init_state ! initialize the ice state + call init_transport ! initialize horizontal transport + call ice_HaloRestore_init ! restored boundary conditions + + call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & + wave_spec_out=wave_spec) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (skl_bgc .or. z_tracers) call alloc_forcing_bgc ! allocate biogeochemistry arrays + + call init_restart ! initialize restart variables + call init_diags ! initialize diagnostic output points + call init_history_therm ! initialize thermo history variables + call init_history_dyn ! initialize dynamic history variables + + call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) + call icepack_query_tracer_flags(tr_iso_out=tr_iso) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (tr_aero .or. tr_zaero) then + call faero_optics !initialize aerosol optical property tables + end if + + ! Initialize shortwave components using swdn from previous timestep + ! if restarting. These components will be scaled to current forcing + ! in prep_radiation. + + if (trim(runtype) == 'continue' .or. restart) then + call init_shortwave ! initialize radiative transfer + end if + + !-------------------------------------------------------------------- + ! coupler communication or forcing data initialization + !-------------------------------------------------------------------- + + if (z_tracers) call get_atm_bgc ! biogeochemistry + + if (runtype == 'initial' .and. .not. restart) then + call init_shortwave ! initialize radiative transfer using current swdn + end if + + call init_flux_atm ! initialize atmosphere fluxes sent to coupler + call init_flux_ocn ! initialize ocean fluxes sent to coupler + + end subroutine cice_init2 + + !======================================================================= + + subroutine init_restart() + + use ice_arrays_column, only: dhsn + use ice_blocks, only: nx_block, ny_block + use ice_calendar, only: time, calendar + use ice_constants, only: c0 + use ice_domain, only: nblocks + use ice_domain_size, only: ncat, n_iso, n_aero, nfsd + use ice_dyn_eap, only: read_restart_eap + use ice_dyn_shared, only: kdyn + use ice_grid, only: tmask + use ice_init, only: ice_ic + use ice_init_column, only: init_age, init_FY, init_lvl, & + init_meltponds_cesm, init_meltponds_lvl, init_meltponds_topo, & + init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd + use ice_restart_column, only: restart_age, read_restart_age, & + restart_FY, read_restart_FY, restart_lvl, read_restart_lvl, & + restart_pond_cesm, read_restart_pond_cesm, & + restart_pond_lvl, read_restart_pond_lvl, & + restart_pond_topo, read_restart_pond_topo, & + restart_fsd, read_restart_fsd, & + restart_iso, read_restart_iso, & + restart_aero, read_restart_aero, & + restart_hbrine, read_restart_hbrine, & + restart_zsal, restart_bgc + use ice_restart_driver, only: restartfile + use ice_restart_shared, only: runtype, restart + use ice_state ! almost everything + + integer(kind=int_kind) :: & i, j , & ! horizontal indices iblk ! block index - logical(kind=log_kind) :: & - tr_iage, tr_FY, tr_lvl, tr_pond_cesm, tr_pond_lvl, & - tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, & - skl_bgc, z_tracers, solve_zsal - integer(kind=int_kind) :: & - ntrcr - integer(kind=int_kind) :: & - nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & - nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice - - character(len=*), parameter :: subname = '(init_restart)' - - call icepack_query_tracer_sizes(ntrcr_out=ntrcr) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & - file=__FILE__, line=__LINE__) - - call icepack_query_parameters(skl_bgc_out=skl_bgc, & - z_tracers_out=z_tracers, solve_zsal_out=solve_zsal) - call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & - tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & - tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & - tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) - call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & - nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & - nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, & - nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + logical(kind=log_kind) :: & + tr_iage, tr_FY, tr_lvl, tr_pond_cesm, tr_pond_lvl, & + tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, & + skl_bgc, z_tracers, solve_zsal + integer(kind=int_kind) :: & + ntrcr + integer(kind=int_kind) :: & + nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & + nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice + + character(len=*), parameter :: subname = '(init_restart)' + !---------------------------------------------------- + + call icepack_query_tracer_sizes(ntrcr_out=ntrcr) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - if (trim(runtype) == 'continue') then - ! start from core restart file - call restartfile() ! given by pointer in ice_in - call calendar(time) ! update time parameters - if (kdyn == 2) call read_restart_eap ! EAP - else if (restart) then ! ice_ic = core restart file - call restartfile (ice_ic) ! or 'default' or 'none' - !!! uncomment to create netcdf - ! call restartfile_v4 (ice_ic) ! CICE v4.1 binary restart file - !!! uncomment if EAP restart data exists - ! if (kdyn == 2) call read_restart_eap - endif - - ! tracers - ! ice age tracer - if (tr_iage) then - if (trim(runtype) == 'continue') & - restart_age = .true. - if (restart_age) then - call read_restart_age - else - do iblk = 1, nblocks - call init_age(trcrn(:,:,nt_iage,:,iblk)) - enddo ! iblk - endif - endif - ! first-year area tracer - if (tr_FY) then - if (trim(runtype) == 'continue') restart_FY = .true. - if (restart_FY) then - call read_restart_FY - else - do iblk = 1, nblocks - call init_FY(trcrn(:,:,nt_FY,:,iblk)) - enddo ! iblk - endif - endif - ! level ice tracer - if (tr_lvl) then - if (trim(runtype) == 'continue') restart_lvl = .true. - if (restart_lvl) then - call read_restart_lvl - else - do iblk = 1, nblocks - call init_lvl(iblk,trcrn(:,:,nt_alvl,:,iblk), & - trcrn(:,:,nt_vlvl,:,iblk)) - enddo ! iblk - endif - endif - ! CESM melt ponds - if (tr_pond_cesm) then - if (trim(runtype) == 'continue') & - restart_pond_cesm = .true. - if (restart_pond_cesm) then - call read_restart_pond_cesm - else - do iblk = 1, nblocks - call init_meltponds_cesm(trcrn(:,:,nt_apnd,:,iblk), & - trcrn(:,:,nt_hpnd,:,iblk)) - enddo ! iblk - endif - endif - ! level-ice melt ponds - if (tr_pond_lvl) then - if (trim(runtype) == 'continue') & - restart_pond_lvl = .true. - if (restart_pond_lvl) then - call read_restart_pond_lvl - else - do iblk = 1, nblocks - call init_meltponds_lvl(trcrn(:,:,nt_apnd,:,iblk), & - trcrn(:,:,nt_hpnd,:,iblk), & - trcrn(:,:,nt_ipnd,:,iblk), & - dhsn(:,:,:,iblk)) - enddo ! iblk - endif - endif - ! topographic melt ponds - if (tr_pond_topo) then - if (trim(runtype) == 'continue') & - restart_pond_topo = .true. - if (restart_pond_topo) then - call read_restart_pond_topo - else - do iblk = 1, nblocks - call init_meltponds_topo(trcrn(:,:,nt_apnd,:,iblk), & - trcrn(:,:,nt_hpnd,:,iblk), & - trcrn(:,:,nt_ipnd,:,iblk)) - enddo ! iblk - endif ! .not. restart_pond - endif - ! floe size distribution - if (tr_fsd) then - if (trim(runtype) == 'continue') restart_fsd = .true. - if (restart_fsd) then - call read_restart_fsd - else - call init_fsd(trcrn(:,:,nt_fsd:nt_fsd+nfsd-1,:,:)) - endif - endif - - ! isotopes - if (tr_iso) then - if (trim(runtype) == 'continue') restart_iso = .true. - if (restart_iso) then - call read_restart_iso - else - do iblk = 1, nblocks - call init_isotope(trcrn(:,:,nt_isosno:nt_isosno+n_iso-1,:,iblk), & - trcrn(:,:,nt_isoice:nt_isoice+n_iso-1,:,iblk)) - enddo ! iblk - endif - endif - - if (tr_aero) then ! ice aerosol - if (trim(runtype) == 'continue') restart_aero = .true. - if (restart_aero) then - call read_restart_aero - else - do iblk = 1, nblocks - call init_aerosol(trcrn(:,:,nt_aero:nt_aero+4*n_aero-1,:,iblk)) - enddo ! iblk - endif ! .not. restart_aero - endif - - if (trim(runtype) == 'continue') then - if (tr_brine) & - restart_hbrine = .true. - if (solve_zsal) & - restart_zsal = .true. - if (skl_bgc .or. z_tracers) & - restart_bgc = .true. - endif - - if (tr_brine .or. skl_bgc) then ! brine height tracer - call init_hbrine - if (tr_brine .and. restart_hbrine) call read_restart_hbrine - endif - - if (solve_zsal .or. skl_bgc .or. z_tracers) then ! biogeochemistry - if (tr_fsd) then - write (nu_diag,*) 'FSD implementation incomplete for use with BGC' - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + call icepack_query_parameters(skl_bgc_out=skl_bgc, & + z_tracers_out=z_tracers, solve_zsal_out=solve_zsal) + call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & + tr_lvl_out=tr_lvl, tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & + tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & + tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) + call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & + nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & + nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, & + nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + if (trim(runtype) == 'continue') then + ! start from core restart file + call restartfile() ! given by pointer in ice_in + call calendar(time) ! update time parameters + if (kdyn == 2) call read_restart_eap ! EAP + else if (restart) then ! ice_ic = core restart file + call restartfile (ice_ic) ! or 'default' or 'none' +!!! uncomment to create netcdf + ! call restartfile_v4 (ice_ic) ! CICE v4.1 binary restart file +!!! uncomment if EAP restart data exists + ! if (kdyn == 2) call read_restart_eap + endif + + ! tracers + ! ice age tracer + if (tr_iage) then + if (trim(runtype) == 'continue') & + restart_age = .true. + if (restart_age) then + call read_restart_age + else + do iblk = 1, nblocks + call init_age(trcrn(:,:,nt_iage,:,iblk)) + enddo ! iblk + endif + endif + ! first-year area tracer + if (tr_FY) then + if (trim(runtype) == 'continue') restart_FY = .true. + if (restart_FY) then + call read_restart_FY + else + do iblk = 1, nblocks + call init_FY(trcrn(:,:,nt_FY,:,iblk)) + enddo ! iblk + endif + endif + ! level ice tracer + if (tr_lvl) then + if (trim(runtype) == 'continue') restart_lvl = .true. + if (restart_lvl) then + call read_restart_lvl + else + do iblk = 1, nblocks + call init_lvl(iblk,trcrn(:,:,nt_alvl,:,iblk), & + trcrn(:,:,nt_vlvl,:,iblk)) + enddo ! iblk + endif + endif + ! CESM melt ponds + if (tr_pond_cesm) then + if (trim(runtype) == 'continue') & + restart_pond_cesm = .true. + if (restart_pond_cesm) then + call read_restart_pond_cesm + else + do iblk = 1, nblocks + call init_meltponds_cesm(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk)) + enddo ! iblk + endif + endif + ! level-ice melt ponds + if (tr_pond_lvl) then + if (trim(runtype) == 'continue') & + restart_pond_lvl = .true. + if (restart_pond_lvl) then + call read_restart_pond_lvl + else + do iblk = 1, nblocks + call init_meltponds_lvl(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk), & + trcrn(:,:,nt_ipnd,:,iblk), & + dhsn(:,:,:,iblk)) + enddo ! iblk + endif + endif + ! topographic melt ponds + if (tr_pond_topo) then + if (trim(runtype) == 'continue') & + restart_pond_topo = .true. + if (restart_pond_topo) then + call read_restart_pond_topo + else + do iblk = 1, nblocks + call init_meltponds_topo(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk), & + trcrn(:,:,nt_ipnd,:,iblk)) + enddo ! iblk + endif ! .not. restart_pond + endif + ! floe size distribution + if (tr_fsd) then + if (trim(runtype) == 'continue') restart_fsd = .true. + if (restart_fsd) then + call read_restart_fsd + else + call init_fsd(trcrn(:,:,nt_fsd:nt_fsd+nfsd-1,:,:)) + endif + endif + + ! isotopes + if (tr_iso) then + if (trim(runtype) == 'continue') restart_iso = .true. + if (restart_iso) then + call read_restart_iso + else + do iblk = 1, nblocks + call init_isotope(trcrn(:,:,nt_isosno:nt_isosno+n_iso-1,:,iblk), & + trcrn(:,:,nt_isoice:nt_isoice+n_iso-1,:,iblk)) + enddo ! iblk + endif + endif + + if (tr_aero) then ! ice aerosol + if (trim(runtype) == 'continue') restart_aero = .true. + if (restart_aero) then + call read_restart_aero + else + do iblk = 1, nblocks + call init_aerosol(trcrn(:,:,nt_aero:nt_aero+4*n_aero-1,:,iblk)) + enddo ! iblk + endif ! .not. restart_aero + endif + + if (trim(runtype) == 'continue') then + if (tr_brine) & + restart_hbrine = .true. + if (solve_zsal) & + restart_zsal = .true. + if (skl_bgc .or. z_tracers) & + restart_bgc = .true. + endif + + if (tr_brine .or. skl_bgc) then ! brine height tracer + call init_hbrine + if (tr_brine .and. restart_hbrine) call read_restart_hbrine + endif + + if (solve_zsal .or. skl_bgc .or. z_tracers) then ! biogeochemistry + if (tr_fsd) then + write (nu_diag,*) 'FSD implementation incomplete for use with BGC' + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - endif - call init_bgc - endif - - !----------------------------------------------------------------- - ! aggregate tracers - !----------------------------------------------------------------- - - !$OMP PARALLEL DO PRIVATE(iblk) - do iblk = 1, nblocks - do j = 1, ny_block - do i = 1, nx_block - if (tmask(i,j,iblk)) then - call icepack_aggregate(ncat = ncat, & - aicen = aicen(i,j,:,iblk), & - trcrn = trcrn(i,j,:,:,iblk), & - vicen = vicen(i,j,:,iblk), & - vsnon = vsnon(i,j,:,iblk), & - aice = aice (i,j, iblk), & - trcr = trcr (i,j,:,iblk), & - vice = vice (i,j, iblk), & - vsno = vsno (i,j, iblk), & - aice0 = aice0(i,j, iblk), & - ntrcr = ntrcr, & - trcr_depend = trcr_depend, & - trcr_base = trcr_base, & - n_trcr_strata = n_trcr_strata, & - nt_strata = nt_strata) - else - ! tcraig, reset all tracer values on land to zero - trcrn(i,j,:,:,iblk) = c0 - endif - enddo - enddo - enddo - !$OMP END PARALLEL DO - - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + endif + call init_bgc + endif + + !----------------------------------------------------------------- + ! aggregate tracers + !----------------------------------------------------------------- + + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + if (tmask(i,j,iblk)) then + call icepack_aggregate(ncat = ncat, & + aicen = aicen(i,j,:,iblk), & + trcrn = trcrn(i,j,:,:,iblk), & + vicen = vicen(i,j,:,iblk), & + vsnon = vsnon(i,j,:,iblk), & + aice = aice (i,j, iblk), & + trcr = trcr (i,j,:,iblk), & + vice = vice (i,j, iblk), & + vsno = vsno (i,j, iblk), & + aice0 = aice0(i,j, iblk), & + ntrcr = ntrcr, & + trcr_depend = trcr_depend, & + trcr_base = trcr_base, & + n_trcr_strata = n_trcr_strata, & + nt_strata = nt_strata) + else + ! tcraig, reset all tracer values on land to zero + trcrn(i,j,:,:,iblk) = c0 + endif + enddo + enddo + enddo + !$OMP END PARALLEL DO + + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - end subroutine init_restart + end subroutine init_restart -!======================================================================= + !======================================================================= - end module CICE_InitMod +end module CICE_InitMod !======================================================================= diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index b70e3417b..74e3307ec 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -83,7 +83,7 @@ module ice_comp_nuopc character(len=*),parameter :: shr_cal_noleap = 'NO_LEAP' character(len=*),parameter :: shr_cal_gregorian = 'GREGORIAN' - type(ESMF_Mesh) :: Emesh + type(ESMF_Mesh) :: ice_mesh integer :: dbug = 0 integer , parameter :: debug_import = 0 ! internal debug level @@ -178,22 +178,16 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Local variables character(len=char_len_long) :: cvalue - character(len=char_len_long) :: logmsg logical :: isPresent, isSet real(dbl_kind) :: eccen, obliqr, lambm0, mvelpp - type(ESMF_DistGrid) :: distGrid - type(ESMF_Mesh) :: EmeshTemp - integer :: spatialDim - integer :: numOwnedElements - real(dbl_kind), pointer :: ownedElemCoords(:) - real(dbl_kind), pointer :: lat(:), latMesh(:) - real(dbl_kind), pointer :: lon(:), lonMesh(:) - integer , allocatable :: gindex_ice(:) - integer , allocatable :: gindex_elim(:) - integer , allocatable :: gindex(:) - integer :: globalID + type(ESMF_DistGrid) :: ice_distGrid character(len=char_len) :: tfrz_option + integer(int_kind) :: ktherm + integer :: localPet + integer :: npes + logical :: mastertask type(ESMF_VM) :: vm + integer :: lmpicom ! local communicator type(ESMF_Time) :: currTime ! Current time type(ESMF_Time) :: startTime ! Start time type(ESMF_Time) :: stopTime ! Stop time @@ -212,31 +206,14 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) integer :: yy,mm,dd ! Temporaries for time query integer :: iyear ! yyyy integer :: dtime ! time step - integer :: lmpicom integer :: shrlogunit ! original log unit character(len=char_len) :: starttype ! infodata start type integer :: lsize ! local size of coupling array - integer :: localPet integer :: n,c,g,i,j,m ! indices integer :: iblk, jblk ! indices integer :: ig, jg ! indices integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain type(block) :: this_block ! block information for current block - integer :: compid ! component id - character(len=char_len_long) :: tempc1,tempc2 - real(dbl_kind) :: diff_lon - integer :: npes - integer :: num_elim_global - integer :: num_elim_local - integer :: num_elim - integer :: num_ice - integer :: num_elim_gcells ! local number of eliminated gridcells - integer :: num_elim_blocks ! local number of eliminated blocks - integer :: num_total_blocks - integer :: my_elim_start, my_elim_end - real(dbl_kind) :: rad_to_deg - integer(int_kind) :: ktherm - logical :: mastertask character(len=char_len_long) :: diag_filename = 'unset' character(len=*), parameter :: subname=trim(modName)//':(InitializeAdvertise) ' !-------------------------------- @@ -515,15 +492,60 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) end if !---------------------------------------------------------------------------- - ! Initialize cice + ! First cice initialization phase - before initializing grid info !---------------------------------------------------------------------------- - ! Note that cice_init also sets time manager info as well as mpi communicator info, - ! including master_task and my_task + call t_startf ('cice_init1') + call cice_init1 + call t_stopf ('cice_init1') + + !---------------------------------------------------------------------------- + ! Initialize grid info + !---------------------------------------------------------------------------- + + ! Determine the model distgrid using the decomposition obtained in + ! call to init_grid1 called from cice_init1 + call ice_mesh_set_distgrid(localpet, npes, ice_distgrid, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! read in the mesh + call NUOPC_CompAttributeGet(gcomp, name='mesh_ice', value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (my_task == master_task) then + write(nu_diag,*)'mesh file for cice domain is ',trim(cvalue) + end if - call t_startf ('cice_init') - call cice_init - call t_stopf ('cice_init') + ! Read in the ice mesh on the cice distribution + ice_mesh = ESMF_MeshCreate(filename=trim(ice_meshfile), fileformat=ESMF_FILEFORMAT_ESMFMESH, & + elementDistGrid=ice_distgrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Initialize the cice mesh and the cice mask + if (trim(grid_type) == 'setmask') then + ! Determine mask input file and create the mask + call NUOPC_CompAttributeGet(gcomp, name='mesh_mask', value=ice_maskfile, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (my_task == master_task) then + write(nu_diag,*)'mask file for cice domain is ',trim(ice_maskfile) + end if + call ice_mesh_setmask_from_maskfile(ice_mesh, ice_maskfile, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + ! In this case init_grid2 will initialize tlon, tlat, area and hm + call init_grid2() + call ice_mesh_check(ice_mesh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + !---------------------------------------------------------------------------- + ! Second cice initialization phase -after initializing grid info + !---------------------------------------------------------------------------- + + ! Note that cice_init2 also sets time manager info as well as mpi communicator info, + ! including master_task and my_task + call t_startf ('cice_init2') + call cice_init2() + call t_stopf ('cice_init2') !---------------------------------------------------------------------------- ! reset shr logging to my log file @@ -611,218 +633,11 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call accum_hist(dt) ! write initial conditions end if - !--------------------------------------------------------------------------- - ! Determine the global index space needed for the distgrid - !--------------------------------------------------------------------------- - - ! number the local grid to get allocation size for gindex_ice - lsize = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - lsize = lsize + 1 - enddo - enddo - enddo - - ! set global index array - allocate(gindex_ice(lsize)) - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) - gindex_ice(n) = (jg-1)*nx_global + ig - enddo - enddo - enddo - - ! Determine total number of eliminated blocks globally - globalID = 0 - num_elim_global = 0 ! number of eliminated blocks - num_total_blocks = 0 - do jblk=1,nblocks_y - do iblk=1,nblocks_x - globalID = globalID + 1 - num_total_blocks = num_total_blocks + 1 - if (distrb_info%blockLocation(globalID) == 0) then - num_elim_global = num_elim_global + 1 - end if - end do - end do - - if (num_elim_global > 0) then - - ! Distribute the eliminated blocks in a round robin fashion amoung processors - num_elim_local = num_elim_global / npes - my_elim_start = num_elim_local*localPet + min(localPet, mod(num_elim_global, npes)) + 1 - if (localPet < mod(num_elim_global, npes)) then - num_elim_local = num_elim_local + 1 - end if - my_elim_end = my_elim_start + num_elim_local - 1 - - ! Determine the number of eliminated gridcells locally - globalID = 0 - num_elim_blocks = 0 ! local number of eliminated blocks - num_elim_gcells = 0 - do jblk=1,nblocks_y - do iblk=1,nblocks_x - globalID = globalID + 1 - if (distrb_info%blockLocation(globalID) == 0) then - num_elim_blocks = num_elim_blocks + 1 - if (num_elim_blocks >= my_elim_start .and. num_elim_blocks <= my_elim_end) then - this_block = get_block(globalID, globalID) - num_elim_gcells = num_elim_gcells + & - (this_block%jhi-this_block%jlo+1) * (this_block%ihi-this_block%ilo+1) - end if - end if - end do - end do - - ! Determine the global index space of the eliminated gridcells - allocate(gindex_elim(num_elim_gcells)) - globalID = 0 - num_elim_gcells = 0 ! local number of eliminated gridcells - num_elim_blocks = 0 ! local number of eliminated blocks - do jblk=1,nblocks_y - do iblk=1,nblocks_x - globalID = globalID + 1 - if (distrb_info%blockLocation(globalID) == 0) then - this_block = get_block(globalID, globalID) - num_elim_blocks = num_elim_blocks + 1 - if (num_elim_blocks >= my_elim_start .and. num_elim_blocks <= my_elim_end) then - do j=this_block%jlo,this_block%jhi - do i=this_block%ilo,this_block%ihi - num_elim_gcells = num_elim_gcells + 1 - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) - gindex_elim(num_elim_gcells) = (jg-1)*nx_global + ig - end do - end do - end if - end if - end do - end do - - ! create a global index that includes both active and eliminated gridcells - num_ice = size(gindex_ice) - num_elim = size(gindex_elim) - allocate(gindex(num_elim + num_ice)) - do n = 1,num_ice - gindex(n) = gindex_ice(n) - end do - do n = num_ice+1,num_ice+num_elim - gindex(n) = gindex_elim(n-num_ice) - end do - - deallocate(gindex_elim) - - else - - ! No eliminated land blocks - num_ice = size(gindex_ice) - allocate(gindex(num_ice)) - do n = 1,num_ice - gindex(n) = gindex_ice(n) - end do - - end if - - !--------------------------------------------------------------------------- - ! Create distGrid from global index array - !--------------------------------------------------------------------------- - - DistGrid = ESMF_DistGridCreate(arbSeqIndexList=gindex, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - !--------------------------------------------------------------------------- - ! Create the CICE mesh - !--------------------------------------------------------------------------- - - ! read in the mesh - call NUOPC_CompAttributeGet(gcomp, name='mesh_ice', value=cvalue, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - EMeshTemp = ESMF_MeshCreate(filename=trim(cvalue), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (my_task == master_task) then - write(nu_diag,*)'mesh file for cice domain is ',trim(cvalue) - end if - - ! recreate the mesh using the above distGrid - EMesh = ESMF_MeshCreate(EMeshTemp, elementDistgrid=Distgrid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - ! obtain mesh lats and lons - call ESMF_MeshGet(Emesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - allocate(ownedElemCoords(spatialDim*numOwnedElements)) - allocate(lonMesh(numOwnedElements), latMesh(numOwnedElements)) - call ESMF_MeshGet(Emesh, ownedElemCoords=ownedElemCoords) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - do n = 1,numOwnedElements - lonMesh(n) = ownedElemCoords(2*n-1) - latMesh(n) = ownedElemCoords(2*n) - end do - - call icepack_query_parameters(rad_to_deg_out=rad_to_deg) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & - file=__FILE__, line=__LINE__) - - ! obtain internally generated cice lats and lons for error checks - allocate(lon(lsize)) - allocate(lat(lsize)) - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - lon(n) = tlon(i,j,iblk)*rad_to_deg - lat(n) = tlat(i,j,iblk)*rad_to_deg - enddo - enddo - enddo - - ! error check differences between internally generated lons and those read in - do n = 1,lsize - diff_lon = abs(lonMesh(n) - lon(n)) - if ( (diff_lon > 1.e2 .and. abs(diff_lon - 360_dbl_kind) > 1.e-1) .or.& - (diff_lon > 1.e-3 .and. diff_lon < 1._dbl_kind) ) then - !write(6,100)n,lonMesh(n),lon(n), diff_lon -100 format('ERROR: CICE n, lonmesh(n), lon(n), diff_lon = ',i6,2(f21.13,3x),d21.5) - !call abort_ice() - end if - if (abs(latMesh(n) - lat(n)) > 1.e-1) then - !write(6,101)n,latMesh(n),lat(n), abs(latMesh(n)-lat(n)) -101 format('ERROR: CICE n, latmesh(n), lat(n), diff_lat = ',i6,2(f21.13,3x),d21.5) - !call abort_ice() - end if - end do - !----------------------------------------------------------------- ! Prescribed ice initialization !----------------------------------------------------------------- - call ice_prescribed_init(clock, Emesh, rc) + call ice_prescribed_init(clock, ice_mesh, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return !----------------------------------------------------------------- @@ -835,16 +650,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call ice_advertise_fields(gcomp, importState, exportState, flds_scalar_name, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - !----------------------------------------------------------------- - ! deallocate memory - !----------------------------------------------------------------- - - deallocate(ownedElemCoords) - deallocate(lon, lonMesh) - deallocate(lat, latMesh) - deallocate(gindex_ice) - deallocate(gindex) - call t_stopf ('cice_init_total') end subroutine InitializeAdvertise @@ -871,7 +676,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! Realize the actively coupled fields !----------------------------------------------------------------- - call ice_realize_fields(gcomp, mesh=Emesh, & + call ice_realize_fields(gcomp, mesh=ice_mesh, & flds_scalar_name=flds_scalar_name, flds_scalar_num=flds_scalar_num, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return From 66356cb3233eb387d1b3fbee6094619228c74c6d Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Thu, 28 Jan 2021 20:13:18 -0700 Subject: [PATCH 2/4] added ice_mesh_mod --- cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 | 704 ++++++++++++++++++ 1 file changed, 704 insertions(+) create mode 100644 cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 diff --git a/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 b/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 new file mode 100644 index 000000000..b02865320 --- /dev/null +++ b/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 @@ -0,0 +1,704 @@ +module ice_mesh_mod + + use ESMF + use ice_kinds_mod , only : dbl_kind, int_kind + use ice_domain_size , only : nx_global, ny_global, max_blocks + use ice_domain , only : nblocks, blocks_ice, distrb_info + use ice_blocks , only : block, get_block, nx_block, ny_block, nblocks_x, nblocks_y + use ice_shr_methods , only : chkerr + use ice_fileunits , only : nu_diag + use ice_communicate , only : my_task, master_task + use ice_initfc , only : icepack_query_parameters + implicit none + private + + public :: ice_mesh_set_distgrid + public :: ice_mesh_setmask_from_maskfile + public :: ice_mesh_check + + private :: ice_mesh_create_mask + + ! Only relevant for lat-lon grids gridcell value of [1 - (land fraction)] (T-cell) + real (dbl_kind), allocatable, public :: ocn_gridcell_frac(:,:,:) + + character(*), parameter :: u_FILE_u = & + __FILE__ + +!======================================================================= +contains +!======================================================================= + + subroutine ice_mesh_set_distgrid(localpet, npes, distgrid, rc) + + ! Determine the global index space needed for the distgrid + + ! input/output variables + integer , intent(in) :: localpet + integer , intent(in) :: npes + type(ESMF_DistGrid) , intent(inout) :: distgrid + integer , intent(out) :: rc + + ! local variables + integer :: n,c,g,i,j,m ! indices + integer :: iblk, jblk ! indices + integer :: ig, jg ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + integer :: lsize ! local size of coupling array + type(block) :: this_block ! block information for current block + integer :: num_elim_global + integer :: num_elim_local + integer :: num_elim + integer :: num_ice + integer :: num_elim_gcells ! local number of eliminated gridcells + integer :: num_elim_blocks ! local number of eliminated blocks + integer :: num_total_blocks + integer :: my_elim_start, my_elim_end + integer , allocatable :: gindex(:) + integer , allocatable :: gindex_ice(:) + integer , allocatable :: gindex_elim(:) + integer :: globalID + character(len=*), parameter :: subname = ' ice_mesh_set_distgrid: ' + !---------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! number the local grid to get allocation size for gindex_ice + lsize = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + lsize = lsize + 1 + enddo + enddo + enddo + + ! set global index array + allocate(gindex_ice(lsize)) + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + ig = this_block%i_glob(i) + jg = this_block%j_glob(j) + gindex_ice(n) = (jg-1)*nx_global + ig + enddo + enddo + enddo + + ! Determine total number of eliminated blocks globally + globalID = 0 + num_elim_global = 0 ! number of eliminated blocks + num_total_blocks = 0 + do jblk=1,nblocks_y + do iblk=1,nblocks_x + globalID = globalID + 1 + num_total_blocks = num_total_blocks + 1 + if (distrb_info%blockLocation(globalID) == 0) then + num_elim_global = num_elim_global + 1 + end if + end do + end do + + if (num_elim_global > 0) then + + ! Distribute the eliminated blocks in a round robin fashion amoung processors + num_elim_local = num_elim_global / npes + my_elim_start = num_elim_local*localPet + min(localPet, mod(num_elim_global, npes)) + 1 + if (localPet < mod(num_elim_global, npes)) then + num_elim_local = num_elim_local + 1 + end if + my_elim_end = my_elim_start + num_elim_local - 1 + + ! Determine the number of eliminated gridcells locally + globalID = 0 + num_elim_blocks = 0 ! local number of eliminated blocks + num_elim_gcells = 0 + do jblk=1,nblocks_y + do iblk=1,nblocks_x + globalID = globalID + 1 + if (distrb_info%blockLocation(globalID) == 0) then + num_elim_blocks = num_elim_blocks + 1 + if (num_elim_blocks >= my_elim_start .and. num_elim_blocks <= my_elim_end) then + this_block = get_block(globalID, globalID) + num_elim_gcells = num_elim_gcells + & + (this_block%jhi-this_block%jlo+1) * (this_block%ihi-this_block%ilo+1) + end if + end if + end do + end do + + ! Determine the global index space of the eliminated gridcells + allocate(gindex_elim(num_elim_gcells)) + globalID = 0 + num_elim_gcells = 0 ! local number of eliminated gridcells + num_elim_blocks = 0 ! local number of eliminated blocks + do jblk=1,nblocks_y + do iblk=1,nblocks_x + globalID = globalID + 1 + if (distrb_info%blockLocation(globalID) == 0) then + this_block = get_block(globalID, globalID) + num_elim_blocks = num_elim_blocks + 1 + if (num_elim_blocks >= my_elim_start .and. num_elim_blocks <= my_elim_end) then + do j=this_block%jlo,this_block%jhi + do i=this_block%ilo,this_block%ihi + num_elim_gcells = num_elim_gcells + 1 + ig = this_block%i_glob(i) + jg = this_block%j_glob(j) + gindex_elim(num_elim_gcells) = (jg-1)*nx_global + ig + end do + end do + end if + end if + end do + end do + + ! create a global index that includes both active and eliminated gridcells + num_ice = size(gindex_ice) + num_elim = size(gindex_elim) + allocate(gindex(num_elim + num_ice)) + do n = 1,num_ice + gindex(n) = gindex_ice(n) + end do + do n = num_ice+1,num_ice+num_elim + gindex(n) = gindex_elim(n-num_ice) + end do + + deallocate(gindex_elim) + + else + + ! No eliminated land blocks + num_ice = size(gindex_ice) + allocate(gindex(num_ice)) + do n = 1,num_ice + gindex(n) = gindex_ice(n) + end do + + end if + + !--------------------------------------------------------------------------- + ! Create distGrid from global index array + !--------------------------------------------------------------------------- + + DistGrid = ESMF_DistGridCreate(arbSeqIndexList=gindex, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + deallocate(gindex_ice) + deallocate(gindex) + + end subroutine ice_mesh_set_distgrid + + !======================================================================= + subroutine ice_mesh_setmask_from_maskfile(ice_mesh, ice_maskfile, rc) + + use ice_scam , only : scmlat, scmlon, single_column + use ice_grid , only : tlon, tlat, hm, tarea, ULON, ULAT, HTN, HTE, ANGLE, ANGLET + use ice_grid , only : uarea, uarear, tarear, tinyarea + use ice_grid , only : dxt, dyt, dxu, dyu, dyhx, dxhy, cyp, cxp, cym, cxm + use ice_grid , only : kmt_file, makemask, tmask + use ice_boundary , only : ice_HaloUpdate + use ice_domain , only : blocks_ice, nblocks, halo_info, distrb_info + use ice_constants , only : c0, c1, c2, p5, radius + use ice_constants , only : field_loc_center, field_type_scalar + use ice_read_write , only : ice_open_nc, ice_close_nc + use ice_exit , only : abort_ice + use netcdf + + ! input/output variables + type(ESMF_Mesh) , intent(in) :: ice_mesh + character(len=*) , intent(in) :: ice_maskfile + integer , intent(out) :: rc + + ! local variables + integer :: i,j,n + integer :: iblk, jblk ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + type(ESMF_Field) :: areaField + real(dbl_kind) , pointer :: mesh_areas(:) + integer :: numownedelements + real(dbl_kind) , pointer :: ownedElemCoords(:) + integer :: spatialDim + integer (int_kind) :: ni, nj, ncid + integer (int_kind) :: dimid, varid, ier + type (block) :: this_block ! block information for current block + real (dbl_kind) :: closelat ! Single-column latitude value + real (dbl_kind) :: closelon ! Single-column longitude value + real (dbl_kind) :: closelatidx ! Single-column latitude index to retrieve + real (dbl_kind) :: closelonidx ! Single-column longitude index to retrieve + integer (int_kind) :: start(2) ! Start index to read in + integer (int_kind) :: count(2) ! Number of points to read in + integer (int_kind) :: start3(3) ! Start index to read in + integer (int_kind) :: count3(3) ! Number of points to read in + integer (int_kind) :: status ! status flag + real (dbl_kind), allocatable :: lats(:) ! temporary + real (dbl_kind), allocatable :: lons(:) ! temporary + real (dbl_kind), allocatable :: pos_lons(:) ! temporary + real (dbl_kind), allocatable :: glob_grid(:,:) ! temporary + real (dbl_kind) :: pos_scmlon ! temporary + real (dbl_kind) :: scamdata ! temporary + integer (int_kind), pointer :: ice_mask(:) + real(dbl_kind) , pointer :: ice_frac(:) + real(dbl_kind) :: pi + real(dbl_kind) :: puny + real(dbl_kind) :: deg_to_rad + character(len=*), parameter :: subname = ' ice_mesh_setmask_from_maskfile' + !--------------------------------------------------- + + rc = ESMF_SUCCESS + + ! Determine start/count to read in for either single column or global lat-lon grid + ! If single_column, then assume that only master_task is used since there is only one task + + if (.not. single_column) then + + ! Obtain the model mask and model frac by mapping the mesh created by reading + ! in the model_maskfile to the model mesh and then resetting the model mesh mask + call ice_mesh_create_mask(ice_mesh, ice_maskfile, ice_mask, ice_frac, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Obtain mesh areas in radians^2 + areaField = ESMF_FieldCreate(ice_mesh, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldRegridGetArea(areaField, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(areaField, farrayPtr=mesh_areas, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Obtain mesh lons and lats in degrees + call ESMF_MeshGet(ice_mesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(ownedElemCoords(spatialDim*numownedelements)) + call ESMF_MeshGet(ice_mesh, ownedElemCoords=ownedElemCoords) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_MeshGet(ice_mesh, ownedElemCoords=ownedElemCoords, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Allocate module variable ocn_gridcell_frac + allocate(ocn_gridcell_frac(nx_block,ny_block,max_blocks)) + + ! Get required constants + call icepack_query_parameters(pi_out=pi, pi2_out=pi2, puny_out=puny, deg_to_rad_out=deg_to_rad) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + ! Set tlon, tlat, tarea, hm and ocn_gridcell_frac + ! Convert mesh areas from radians^2 to m^2 (tarea is in m^2) + ! Convert lons and lats from degrees to radians + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n + 1 + tlon(i,j,iblk) = ownedElemCoords(2*n-1) * deg_to_rad + tlat(i,j,iblk) = ownedElemCoords(2*n) * deg_to_rad + tarea(i,j,iblk) = mesh_areas(n) * (radius*radius) + hm(i,j,iblk) = real(ice_mask(n),kind=dbl_kind) + ocn_gridcell_frac(i,j,iblk) = ice_frac(n) + enddo + enddo + enddo + + ! Dealocate memory + deallocate(ownedElemCoords) + call ESMF_FieldDestroy(areaField) + + else ! single column mode + + if (my_task == master_task) then + call ice_open_nc(kmt_file, ncid) + status = nf90_inq_dimid (ncid, 'ni', dimid) + status = nf90_inquire_dimension(ncid, dimid, len=ni) + status = nf90_inq_dimid (ncid, 'nj', dimid) + status = nf90_inquire_dimension(ncid, dimid, len=nj) + end if + + ! Check for consistency + if (my_task == master_task) then + if ((nx_global /= 1).or. (ny_global /= 1)) then + write(nu_diag,*) 'Because you have selected the column model flag' + write(nu_diag,*) 'Please set nx_global=ny_global=1 in file' + write(nu_diag,*) 'ice_domain_size.F and recompile' + call abort_ice ('ice_mesh_setmask_from_maskfile: check nx_global, ny_global') + endif + end if + + ! Read in domain file for single column + allocate(lats(nj)) + allocate(lons(ni)) + allocate(pos_lons(ni)) + allocate(glob_grid(ni,nj)) + + start3=(/1,1,1/) + count3=(/ni,nj,1/) + status = nf90_inq_varid(ncid, 'xc' , varid) + if (status /= nf90_noerr) call abort_ice (subname//' inq_varid xc') + status = nf90_get_var(ncid, varid, glob_grid, start3, count3) + if (status /= nf90_noerr) call abort_ice (subname//' get_var xc') + do i = 1,ni + lons(i) = glob_grid(i,1) + end do + + status = nf90_inq_varid(ncid, 'yc' , varid) + if (status /= nf90_noerr) call abort_ice (subname//' inq_varid yc') + status = nf90_get_var(ncid, varid, glob_grid, start3, count3) + if (status /= nf90_noerr) call abort_ice (subname//' get_var yc') + do j = 1,nj + lats(j) = glob_grid(1,j) + end do + + ! convert lons array and scmlon to 0,360 and find index of value closest to 0 + ! and obtain single-column longitude/latitude indices to retrieve + + pos_lons(:)= mod(lons(:) + 360._dbl_kind,360._dbl_kind) + pos_scmlon = mod(scmlon + 360._dbl_kind,360._dbl_kind) + start(1) = (MINLOC(abs(pos_lons-pos_scmlon),dim=1)) + start(2) = (MINLOC(abs(lats -scmlat ),dim=1)) + + deallocate(lats) + deallocate(lons) + deallocate(pos_lons) + deallocate(glob_grid) + + status = nf90_inq_varid(ncid, 'xc' , varid) + if (status /= nf90_noerr) call abort_ice (subname//' inq_varid xc') + status = nf90_get_var(ncid, varid, scamdata, start) + if (status /= nf90_noerr) call abort_ice (subname//' get_var xc') + TLON = scamdata + status = nf90_inq_varid(ncid, 'yc' , varid) + if (status /= nf90_noerr) call abort_ice (subname//' inq_varid yc') + status = nf90_get_var(ncid, varid, scamdata, start) + if (status /= nf90_noerr) call abort_ice (subname//' get_var yc') + TLAT = scamdata + status = nf90_inq_varid(ncid, 'area' , varid) + if (status /= nf90_noerr) call abort_ice (subname//' inq_varid area') + status = nf90_get_var(ncid, varid, scamdata, start) + if (status /= nf90_noerr) call abort_ice (subname//' get_var are') + tarea = scamdata + status = nf90_inq_varid(ncid, 'mask' , varid) + if (status /= nf90_noerr) call abort_ice (subname//' inq_varid mask') + status = nf90_get_var(ncid, varid, scamdata, start) + if (status /= nf90_noerr) call abort_ice (subname//' get_var mask') + hm = scamdata + status = nf90_inq_varid(ncid, 'frac' , varid) + if (status /= nf90_noerr) call abort_ice (subname//' inq_varid frac') + status = nf90_get_var(ncid, varid, scamdata, start) + if (status /= nf90_noerr) call abort_ice (subname//' get_var frac') + ocn_gridcell_frac = scamdata + + if (my_task == master_task) then + call ice_close_nc(ncid) + end if + + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + + do j = jlo, jhi + do i = ilo, ihi + ! Convert from degrees to radians + TLON(i,j,iblk) = pi*TLON(i,j,iblk)/180._dbl_kind + + ! Convert from degrees to radians + TLAT(i,j,iblk) = pi*TLAT(i,j,iblk)/180._dbl_kind + + ! Convert from radians^2 to m^2 + ! (area in domain file is in radians^2 and tarea is in m^2) + tarea(i,j,iblk) = tarea(i,j,iblk) * (radius*radius) + end do + end do + end do + + end if + + call ice_HaloUpdate (TLON , halo_info, field_loc_center, field_type_scalar, fillValue=c1) + call ice_HaloUpdate (TLAT , halo_info, field_loc_center, field_type_scalar, fillValue=c1) + call ice_HaloUpdate (tarea , halo_info, field_loc_center, field_type_scalar, fillValue=c1) + call ice_HaloUpdate (hm , halo_info, field_loc_center, field_type_scalar, fillValue=c1) + + !----------------------------------------------------------------- + ! CALCULATE various geometric 2d arrays + ! The U grid (velocity) is not used when run with sequential CAM + ! because we only use thermodynamic sea ice. However, ULAT is used + ! in the default initialization of CICE so we calculate it here as + ! a "dummy" so that CICE will initialize with ice. If a no ice + ! initialization is OK (or desired) this can be commented out and + ! ULAT will remain 0 as specified above. ULAT is located at the + ! NE corner of the grid cell, TLAT at the center, so here ULAT is + ! hacked by adding half the latitudinal spacing (in radians) to TLAT. + !----------------------------------------------------------------- + + ANGLET(:,:,:) = c0 + + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + + do j = jlo, jhi + do i = ilo, ihi + + if (ny_global == 1) then + uarea(i,j,iblk) = tarea(i,j, iblk) + else + uarea(i,j,iblk) = p25* & + (tarea(i,j, iblk) + tarea(i+1,j, iblk) & + + tarea(i,j+1,iblk) + tarea(i+1,j+1,iblk)) + endif + tarear(i,j,iblk) = c1/tarea(i,j,iblk) + uarear(i,j,iblk) = c1/uarea(i,j,iblk) + tinyarea(i,j,iblk) = puny*tarea(i,j,iblk) + + if (single_column) then + ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/nj) + else + if (ny_global == 1) then + ULAT (i,j,iblk) = TLAT(i,j,iblk) + else + ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/ny_global) + endif + endif + ULON (i,j,iblk) = c0 + ANGLE (i,j,iblk) = c0 + + HTN (i,j,iblk) = 1.e36_dbl_kind + HTE (i,j,iblk) = 1.e36_dbl_kind + dxt (i,j,iblk) = 1.e36_dbl_kind + dyt (i,j,iblk) = 1.e36_dbl_kind + dxu (i,j,iblk) = 1.e36_dbl_kind + dyu (i,j,iblk) = 1.e36_dbl_kind + dxhy (i,j,iblk) = 1.e36_dbl_kind + dyhx (i,j,iblk) = 1.e36_dbl_kind + cyp (i,j,iblk) = 1.e36_dbl_kind + cxp (i,j,iblk) = 1.e36_dbl_kind + cym (i,j,iblk) = 1.e36_dbl_kind + cxm (i,j,iblk) = 1.e36_dbl_kind + enddo + enddo + enddo + + call ice_HaloUpdate (ULAT, halo_info, field_loc_center, field_type_scalar, fillValue=c1) + + ! Set the boundary values for the T cell land mask (hm) and + ! make the logical land masks for T and U cells (tmask, umask). + ! Also create hemisphere masks (mask-n northern, mask-s southern) + call makemask() + + end subroutine ice_mesh_setmask_from_maskfile + + !=============================================================================== + subroutine ice_mesh_create_mask(ice_mesh, ice_maskfile, ice_mask, ice_frac, rc) + + use ice_constants, only : c0, c1 + + ! input/out variables + type(ESMF_Mesh) , intent(in) :: ice_mesh + character(len=*) , intent(in) :: ice_maskfile + integer , pointer , intent(out) :: ice_mask(:) + real(dbl_kind) , pointer , intent(out) :: ice_frac(:) + integer , intent(out) :: rc + + ! local variables: + type(ESMF_Mesh) :: mesh_mask + type(ESMF_Field) :: field_mask + type(ESMF_Field) :: field_dst + type(ESMF_RouteHandle) :: rhandle + integer :: srcMaskValue = 0 + integer :: dstMaskValue = -987987 ! spval for RH mask values + integer :: srcTermProcessing_Value = 0 + logical :: checkflag = .false. + real(dbl_kind) , pointer :: mask_src(:) ! on mesh created from ice_maskfile + real(dbl_kind) , pointer :: dataptr1d(:) + type(ESMF_DistGrid) :: distgrid_mask + type(ESMF_Array) :: elemMaskArray + integer :: lsize_mask, lsize_dst + integer :: n, spatialDim + real(dbl_kind) :: fminval = 0.001_dbl_kind ! TODO: make this a share constant + real(dbl_kind) :: fmaxval = 1._dbl_kind + real(dbl_kind) :: lfrac + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + mesh_mask = ESMF_MeshCreate(trim(ice_maskfile), fileformat=ESMF_FILEFORMAT_ESMFMESH, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_MeshGet(ice_mesh, spatialDim=spatialDim, numOwnedElements=lsize_dst, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(ice_mask(lsize_dst)) + allocate(ice_frac(lsize_dst)) + + ! create fields on source and destination meshes + field_mask = ESMF_FieldCreate(mesh_mask, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + field_dst = ESMF_FieldCreate(ice_mesh, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! create route handle to map source mask (assume ocean) to destination mesh (assume atm/lnd) + call ESMF_FieldRegridStore(field_mask, field_dst, routehandle=rhandle, & + srcMaskValues=(/srcMaskValue/), dstMaskValues=(/dstMaskValue/), & + regridmethod=ESMF_REGRIDMETHOD_CONSERVE, normType=ESMF_NORMTYPE_DSTAREA, & + srcTermProcessing=srcTermProcessing_Value, & + ignoreDegenerate=.true., unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! fill in values for field_mask with mask on source mesh + call ESMF_MeshGet(mesh_mask, elementdistGrid=distgrid_mask, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_DistGridGet(distgrid_mask, localDe=0, elementCount=lsize_mask, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(mask_src(lsize_mask)) + elemMaskArray = ESMF_ArrayCreate(distgrid_mask, mask_src, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! The following call fills in the values of mask_src + call ESMF_MeshGet(mesh_mask, elemMaskArray=elemMaskArray, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! The following call fills in the values of field_mask + call ESMF_FieldGet(field_mask, farrayptr=dataptr1d, rc=rc) + dataptr1d(:) = mask_src(:) + + ! map source mask to destination mesh - to obtain destination mask and frac + call ESMF_FieldRegrid(field_mask, field_dst, routehandle=rhandle, & + termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call ESMF_MeshGet(ice_mesh, spatialDim=spatialDim, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(field_dst, farrayptr=dataptr1d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! now determine ice_mask and ice_frac + do n = 1,lsize_dst + lfrac = c1 - dataptr1d(n) + if (lfrac > fmaxval) lfrac = c1 + if (lfrac < fminval) lfrac = c0 + ice_frac(n) = c1 - lfrac + if (ice_frac(n) == c0) then + ice_mask(n) = 0 + else + ice_mask(n) = 1 + end if + enddo + + ! reset the model mesh mask + call ESMF_MeshSet(ice_mesh, elementMask=ice_mask, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! deallocate memory + call ESMF_RouteHandleDestroy(rhandle, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldDestroy(field_mask, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldDestroy(field_dst, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + deallocate(mask_src) + + end subroutine ice_mesh_create_mask + + !=============================================================================== + subroutine ice_mesh_check(ice_mesh, rc) + + ! Check CICE mesh + + use ice_constants, only : c1 + use ice_grid , only : tlon, tlat + + ! input/output parameters + type(ESMF_Mesh) , intent(inout) :: ice_mesh + integer , intent(out) :: rc + + ! local variables + type(ESMF_DistGrid) :: distGrid + integer :: n,c,g,i,j,m ! indices + integer :: iblk, jblk ! indices + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + type(block) :: this_block ! block information for current block + integer :: spatialDim + integer :: numOwnedElements + real(dbl_kind), pointer :: ownedElemCoords(:) + real(dbl_kind), pointer :: lat(:), latMesh(:) + real(dbl_kind), pointer :: lon(:), lonMesh(:) + real(dbl_kind) :: diff_lon + real(dbl_kind) :: diff_lat + real(dbl_kind) :: rad_to_deg + !--------------------------------------------------- + + ! error check differences between internally generated lons and those read in + call ESMF_MeshGet(ice_mesh, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(ownedElemCoords(spatialDim*numownedelements)) + allocate(lonmesh(numOwnedElements)) + allocate(latmesh(numOwnedElements)) + call ESMF_MeshGet(ice_mesh, ownedElemCoords=ownedElemCoords) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + do n = 1,numOwnedElements + lonMesh(n) = ownedElemCoords(2*n-1) + latMesh(n) = ownedElemCoords(2*n) + end do + + ! obtain internally generated cice lats and lons for error checks + call icepack_query_parameters(rad_to_deg_out=rad_to_deg) + allocate(lon(numOwnedElements)) + allocate(lat(numOwnedElements)) + lon(:) = 0. + lat(:) = 0. + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n + 1 + lon(n) = tlon(i,j,iblk)*rad_to_deg + lat(n) = tlat(i,j,iblk)*rad_to_deg + + ! error check differences between internally generated lons and those read in + diff_lon = abs(lonMesh(n) - lon(n)) + if ( (diff_lon > 1.e2 .and. abs(diff_lon - 360.) > 1.e-1) .or.& + (diff_lon > 1.e-3 .and. diff_lon < c1) ) then + write(6,100)n,lonMesh(n),lon(n), diff_lon + !call shr_sys_abort() + end if + if (abs(latMesh(n) - lat(n)) > 1.e-1) then + write(6,101)n,latMesh(n),lat(n), abs(latMesh(n)-lat(n)) + !call shr_sys_abort() + end if + + enddo + enddo + enddo + +100 format('ERROR: CICE n, lonmesh, lon, diff_lon = ',i6,2(f21.13,3x),d21.5) +101 format('ERROR: CICE n, latmesh, lat, diff_lat = ',i6,2(f21.13,3x),d21.5) + + ! deallocate memory + deallocate(ownedElemCoords) + deallocate(lon, lonMesh) + deallocate(lat, latMesh) + + end subroutine ice_mesh_check + +end module ice_mesh_mod From 6ff987ca2cd9db47c86a0b9f59a7fe9f54e6299d Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Sun, 31 Jan 2021 18:11:22 -0700 Subject: [PATCH 3/4] implemented area correction factors as option --- cicecore/cicedynB/infrastructure/ice_grid.F90 | 2 +- cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 | 7 +- .../drivers/nuopc/cmeps/ice_comp_nuopc.F90 | 68 +- .../drivers/nuopc/cmeps/ice_import_export.F90 | 919 ++++++++---------- cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 | 57 +- 5 files changed, 476 insertions(+), 577 deletions(-) diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index ce1af48d8..5ebb502cd 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -40,7 +40,7 @@ module ice_grid private public :: init_grid1, init_grid2, & t2ugrid_vector, u2tgrid_vector, & - to_ugrid, to_tgrid, alloc_grid + to_ugrid, to_tgrid, alloc_grid, makemask character (len=char_len_long), public :: & grid_format , & ! file format ('bin'=binary or 'nc'=netcdf) diff --git a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 index a585dd331..7b7c537ce 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 @@ -3,7 +3,7 @@ module CICE_InitMod ! Initialize CICE model. use ice_kinds_mod - use ice_exit, only: abort_ice + use ice_exit , only: abort_ice use ice_fileunits, only: init_fileunits, nu_diag use icepack_intfc, only: icepack_aggregate use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist @@ -15,7 +15,10 @@ module CICE_InitMod implicit none private - public :: cice_init + public :: cice_init1 + public :: cice_init2 + + private :: init_restart !======================================================================= contains diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index 74e3307ec..256fe5bf8 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -20,22 +20,16 @@ module ice_comp_nuopc use ice_shr_methods , only : set_component_logging, get_component_instance, state_flddebug use ice_import_export , only : ice_import, ice_export, ice_advertise_fields, ice_realize_fields use ice_domain_size , only : nx_global, ny_global - use ice_domain , only : nblocks, blocks_ice, distrb_info - use ice_blocks , only : block, get_block, nx_block, ny_block, nblocks_x, nblocks_y - use ice_grid , only : tlon, tlat + use ice_grid , only : grid_type, init_grid2 use ice_communicate , only : init_communicate, my_task, master_task, mpi_comm_ice use ice_calendar , only : force_restart_now, write_ic use ice_calendar , only : idate, mday, time, month, daycal, time2sec, year_init use ice_calendar , only : sec, dt, calendar, calendar_type, nextsw_cday, istep use ice_kinds_mod , only : dbl_kind, int_kind, char_len, char_len_long - use ice_scam , only : scmlat, scmlon, single_column use ice_fileunits , only : nu_diag, nu_diag_set, inst_index, inst_name use ice_fileunits , only : inst_suffix, release_all_fileunits, flush_fileunit use ice_restart_shared , only : runid, runtype, restart, use_restart_time, restart_dir, restart_file use ice_history , only : accum_hist - use ice_flux , only : send_i2x_per_cat - use CICE_InitMod , only : cice_init - use CICE_RunMod , only : cice_run use ice_exit , only : abort_ice use icepack_intfc , only : icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc , only : icepack_init_orbit, icepack_init_parameters, icepack_query_orbit @@ -47,6 +41,9 @@ module ice_comp_nuopc use shr_orb_mod , only : shr_orb_decl, shr_orb_params, SHR_ORB_UNDEF_REAL, SHR_ORB_UNDEF_INT #endif use ice_timers + use CICE_InitMod , only : cice_init1, cice_init2 + use CICE_RunMod , only : cice_run + use ice_mesh_mod , only : ice_mesh_set_distgrid, ice_mesh_setmask_from_maskfile, ice_mesh_check use ice_prescribed_mod , only : ice_prescribed_init implicit none @@ -178,6 +175,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! Local variables character(len=char_len_long) :: cvalue + character(len=char_len_long) :: ice_meshfile logical :: isPresent, isSet real(dbl_kind) :: eccen, obliqr, lambm0, mvelpp type(ESMF_DistGrid) :: ice_distGrid @@ -213,8 +211,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) integer :: iblk, jblk ! indices integer :: ig, jg ! indices integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain - type(block) :: this_block ! block information for current block character(len=char_len_long) :: diag_filename = 'unset' + character(len=char_len_long) :: logmsg character(len=*), parameter :: subname=trim(modName)//':(InitializeAdvertise) ' !-------------------------------- @@ -278,14 +276,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) write(logmsg,'(i6)') dbug call ESMF_LogWrite('CICE_cap: dbug = '//trim(logmsg), ESMF_LOGMSG_INFO) - send_i2x_per_cat = .false. - call NUOPC_CompAttributeGet(gcomp, name='flds_i2o_per_cat', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (isPresent .and. isSet) then - read(cvalue,*) send_i2x_per_cat - call ESMF_LogWrite('send_i2x_per_cat = '// trim(cvalue), ESMF_LOGMSG_INFO) - end if - !---------------------------------------------------------------------------- ! generate local mpi comm !---------------------------------------------------------------------------- @@ -403,23 +393,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) runtype = 'initial' ! determined from the namelist in ice_init if CESMCOUPLED is not defined end if - ! Determine if single column - call NUOPC_CompAttributeGet(gcomp, name='single_column', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (isPresent .and. isSet) then - read(cvalue,*) single_column - if (single_column) then - call NUOPC_CompAttributeGet(gcomp, name='scmlon', value=cvalue, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - read(cvalue,*) scmlon - call NUOPC_CompAttributeGet(gcomp, name='scmlat', value=cvalue, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - read(cvalue,*) scmlat - end if - else - single_column = .false. - end if - ! Determine runid call NUOPC_CompAttributeGet(gcomp, name='case_name', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) if (isPresent .and. isSet) then @@ -474,12 +447,14 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call shr_file_setLogUnit (shrlogunit) - call NUOPC_CompAttributeGet(gcomp, name="diro", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + call NUOPC_CompAttributeGet(gcomp, name="diro", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then diag_filename = trim(cvalue) end if - call NUOPC_CompAttributeGet(gcomp, name="logfile", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + call NUOPC_CompAttributeGet(gcomp, name="logfile", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then diag_filename = trim(diag_filename) // '/' // trim(cvalue) @@ -508,27 +483,19 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call ice_mesh_set_distgrid(localpet, npes, ice_distgrid, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! read in the mesh - call NUOPC_CompAttributeGet(gcomp, name='mesh_ice', value=cvalue, rc=rc) + ! Read in the ice mesh on the cice distribution + call NUOPC_CompAttributeGet(gcomp, name='mesh_ice', value=ice_meshfile, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (my_task == master_task) then - write(nu_diag,*)'mesh file for cice domain is ',trim(cvalue) + write(nu_diag,*)'mesh file for cice domain is ',trim(ice_meshfile) end if - - ! Read in the ice mesh on the cice distribution - ice_mesh = ESMF_MeshCreate(filename=trim(ice_meshfile), fileformat=ESMF_FILEFORMAT_ESMFMESH, & - elementDistGrid=ice_distgrid, rc=rc) + ice_mesh = ESMF_MeshCreate(filename=trim(ice_meshfile), & + fileformat=ESMF_FILEFORMAT_ESMFMESH, elementDistGrid=ice_distgrid, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Initialize the cice mesh and the cice mask if (trim(grid_type) == 'setmask') then - ! Determine mask input file and create the mask - call NUOPC_CompAttributeGet(gcomp, name='mesh_mask', value=ice_maskfile, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (my_task == master_task) then - write(nu_diag,*)'mask file for cice domain is ',trim(ice_maskfile) - end if - call ice_mesh_setmask_from_maskfile(ice_mesh, ice_maskfile, rc=rc) + call ice_mesh_setmask_from_maskfile(gcomp, ice_mesh, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return else ! In this case init_grid2 will initialize tlon, tlat, area and hm @@ -1280,7 +1247,4 @@ subroutine ice_cal_ymd2date(year, month, day, date) end subroutine ice_cal_ymd2date - !=============================================================================== - - end module ice_comp_nuopc diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index a1537a1d8..6e584eb90 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -4,7 +4,7 @@ module ice_import_export use NUOPC use NUOPC_Model use ice_kinds_mod , only : int_kind, dbl_kind, char_len, log_kind - use ice_constants , only : c0, c1, spval_dbl + use ice_constants , only : c0, c1, spval_dbl, radius use ice_constants , only : field_loc_center, field_type_scalar, field_type_vector use ice_blocks , only : block, get_block, nx_block, ny_block use ice_domain , only : nblocks, blocks_ice, halo_info, distrb_info @@ -21,6 +21,7 @@ module ice_import_export use ice_flux , only : fresh, fsalt, zlvl, uatm, vatm, potT, Tair, Qa use ice_flux , only : rhoa, swvdr, swvdf, swidr, swidf, flw, frain use ice_flux , only : fsnow, uocn, vocn, sst, ss_tltx, ss_tlty, frzmlt + use ice_flux , only : send_i2x_per_cat use ice_flux , only : sss, Tf, wind, fsw use ice_state , only : vice, vsno, aice, aicen_init, trcr use ice_grid , only : tlon, tlat, tarea, tmask, anglet, hm, ocn_gridcell_frac @@ -34,7 +35,7 @@ module ice_import_export use icepack_intfc , only : icepack_query_parameters, icepack_query_tracer_flags use icepack_intfc , only : icepack_liquidus_temperature use icepack_intfc , only : icepack_sea_freezing_temperature - use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf + use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf #ifdef CESMCOUPLED use shr_frz_mod , only : shr_frz_freezetemp #endif @@ -54,8 +55,6 @@ module ice_import_export interface state_getfldptr module procedure state_getfldptr_1d module procedure state_getfldptr_2d - module procedure state_getfldptr_3d - module procedure state_getfldptr_4d end interface state_getfldptr private :: state_getfldptr @@ -79,12 +78,15 @@ module ice_import_export integer :: ungridded_ubound = 0 end type fld_list_type + ! area correction factors for fluxes send and received from mediator + real(dbl_kind), allocatable :: mod2med_areacor(:) ! ratios of model areas to input mesh areas + real(dbl_kind), allocatable :: med2mod_areacor(:) ! ratios of input mesh areas to model areas + integer, parameter :: fldsMax = 100 integer :: fldsToIce_num = 0 integer :: fldsFrIce_num = 0 type (fld_list_type) :: fldsToIce(fldsMax) type (fld_list_type) :: fldsFrIce(fldsMax) - type(ESMF_GeomType_Flag) :: geomtype integer , parameter :: io_dbug = 10 ! i/o debug messages character(*), parameter :: u_FILE_u = & @@ -115,13 +117,28 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam rc = ESMF_SUCCESS if (io_dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) + ! Determine if ice sends multiple ice category info back to mediator + send_i2x_per_cat = .false. + call NUOPC_CompAttributeGet(gcomp, name='flds_i2o_per_cat', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) send_i2x_per_cat + end if + if (my_task == master_task) then + write(nu_diag,*)'send_i2x_per_cat = ',send_i2x_per_cat + end if + ! Determine if the following attributes are sent by the driver and if so read them in flds_wiso = .false. - call NUOPC_CompAttributeGet(gcomp, name='flds_wiso', value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc) + call NUOPC_CompAttributeGet(gcomp, name='flds_wiso', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then read(cvalue,*) flds_wiso - call ESMF_LogWrite('flds_wiso = '// trim(cvalue), ESMF_LOGMSG_INFO) + end if + if (my_task == master_task) then + write(nu_diag,*)'flds_wiso = ',flds_wiso end if !----------------- @@ -253,21 +270,27 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam end subroutine ice_advertise_fields -!============================================================================== - - subroutine ice_realize_fields(gcomp, mesh, grid, flds_scalar_name, flds_scalar_num, rc) + !============================================================================== + subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc) ! input/output variables - type(ESMF_GridComp) :: gcomp - type(ESMF_Mesh) , optional , intent(in) :: mesh - type(ESMF_Grid) , optional , intent(in) :: grid - character(len=*) , intent(in) :: flds_scalar_name - integer , intent(in) :: flds_scalar_num - integer , intent(out) :: rc + type(ESMF_GridComp) :: gcomp + type(ESMF_Mesh) , intent(in) :: mesh + character(len=*) , intent(in) :: flds_scalar_name + integer , intent(in) :: flds_scalar_num + integer , intent(out) :: rc ! local variables - type(ESMF_State) :: importState - type(ESMF_State) :: exportState + type(ESMF_State) :: importState + type(ESMF_State) :: exportState + type(ESMF_Field) :: lfield + integer :: numOwnedElements + integer :: i, j, iblk, n + integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + type(block) :: this_block ! block information for current block + real(dbl_kind), allocatable :: mesh_areas(:) + real(dbl_kind), allocatable :: model_areas(:) + real(dbl_kind), pointer :: dataptr(:) character(len=*), parameter :: subname='(ice_import_export:realize_fields)' !--------------------------------------------------------------------------- @@ -276,60 +299,84 @@ subroutine ice_realize_fields(gcomp, mesh, grid, flds_scalar_name, flds_scalar_n call NUOPC_ModelGet(gcomp, importState=importState, exportState=exportState, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - if (present(mesh)) then - - geomtype = ESMF_GEOMTYPE_MESH - - call fldlist_realize( & - state=ExportState, & - fldList=fldsFrIce, & - numflds=fldsFrIce_num, & - flds_scalar_name=flds_scalar_name, & - flds_scalar_num=flds_scalar_num, & - tag=subname//':CICE_Export',& - mesh=mesh, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call fldlist_realize( & - state=importState, & - fldList=fldsToIce, & - numflds=fldsToIce_num, & - flds_scalar_name=flds_scalar_name, & - flds_scalar_num=flds_scalar_num, & - tag=subname//':CICE_Import',& - mesh=mesh, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - else if (present(grid)) then + call fldlist_realize( & + state=ExportState, & + fldList=fldsFrIce, & + numflds=fldsFrIce_num, & + flds_scalar_name=flds_scalar_name, & + flds_scalar_num=flds_scalar_num, & + tag=subname//':CICE_Export',& + mesh=mesh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - geomtype = ESMF_GEOMTYPE_GRID + call fldlist_realize( & + state=importState, & + fldList=fldsToIce, & + numflds=fldsToIce_num, & + flds_scalar_name=flds_scalar_name, & + flds_scalar_num=flds_scalar_num, & + tag=subname//':CICE_Import',& + mesh=mesh, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - call fldlist_realize( & - state=ExportState, & - fldList=fldsFrIce, & - numflds=fldsFrIce_num, & - flds_scalar_name=flds_scalar_name, & - flds_scalar_num=flds_scalar_num, & - tag=subname//':CICE_Export',& - grid=grid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! Determine areas for regridding + call ESMF_MeshGet(mesh, numOwnedElements=numOwnedElements, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate (mod2med_areacor(numOwnedElements)) + allocate (med2mod_areacor(numOwnedElements)) + mod2med_areacor(:) = c1 + med2mod_areacor(:) = c1 - call fldlist_realize( & - state=importState, & - fldList=fldsToIce, & - numflds=fldsToIce_num, & - flds_scalar_name=flds_scalar_name, & - flds_scalar_num=flds_scalar_num, & - tag=subname//':CICE_Import',& - grid=grid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return +#ifdef CESMCOUPLED + ! Get mesh areas from second field - using second field since the + ! first field is the scalar field + call ESMF_StateGet(exportState, itemName=trim(fldsFrIce(2)%stdname), field=lfield, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldRegridGetArea(lfield, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(lfield, farrayPtr=dataptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + allocate(mesh_areas(numOwnedElements)) + mesh_areas(:) = dataptr(:) + + ! Determine model areas + allocate(model_areas(numOwnedElements)) + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + model_areas(n) = tarea(i,j,iblk)/(radius*radius) + enddo + enddo + enddo - end if + ! Determine flux correction factors (module variables) + do n = 1,numOwnedElements + if (model_areas(n) == mesh_areas(n)) then + mod2med_areacor(n) = c1 + med2mod_areacor(n) = c1 + else + mod2med_areacor(n) = model_areas(n) / mesh_areas(n) + med2mod_areacor(n) = mesh_areas(n) / model_areas(n) + if (abs(mod2med_areacor(n) - 1._dbl_kind) > 1.e-13) then + write(6,'(a,i8,2x,d21.14,2x)')' AREACOR cice6: n, abs(mod2med_areacor(n)-1)', & + n, abs(mod2med_areacor(n) - 1._dbl_kind) + end if + end if + end do + deallocate(model_areas) + deallocate(mesh_areas) +#endif end subroutine ice_realize_fields !============================================================================== - subroutine ice_import( importState, rc ) ! input/output variables @@ -346,7 +393,11 @@ subroutine ice_import( importState, rc ) real (kind=dbl_kind) :: workx, worky real (kind=dbl_kind) :: MIN_RAIN_TEMP, MAX_SNOW_TEMP real (kind=dbl_kind) :: Tffresh - real (kind=dbl_kind) :: inst_pres_height_lowest + real (kind=dbl_kind) :: inst_pres_height_lowest + real (kind=dbl_kind), pointer :: dataptr2d(:,:) + real (kind=dbl_kind), pointer :: dataptr1d(:) + real (kind=dbl_kind), pointer :: dataptr2d_dstwet(:,:) + real (kind=dbl_kind), pointer :: dataptr2d_dstdry(:,:) character(len=char_len) :: tfrz_option integer(int_kind) :: ktherm character(len=*), parameter :: subname = 'ice_import' @@ -356,17 +407,20 @@ subroutine ice_import( importState, rc ) call icepack_query_parameters(Tffresh_out=Tffresh) call icepack_query_parameters(tfrz_option_out=tfrz_option) call icepack_query_parameters(ktherm_out=ktherm) - if (io_dbug > 5) then + + if (io_dbug > 5) then write(msgString,'(A,i8)')trim(subname)//' tfrz_option = ' & // trim(tfrz_option)//', ktherm = ',ktherm call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO) end if -! call icepack_query_parameters(tfrz_option_out=tfrz_option, & -! modal_aero_out=modal_aero, z_tracers_out=z_tracers, skl_bgc_out=skl_bgc, & -! Tffresh_out=Tffresh) -! call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_iage_out=tr_iage, & -! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & -! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + + ! call icepack_query_parameters(tfrz_option_out=tfrz_option, & + ! modal_aero_out=modal_aero, z_tracers_out=z_tracers, skl_bgc_out=skl_bgc, & + ! Tffresh_out=Tffresh) + ! call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_iage_out=tr_iage, & + ! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & + ! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=u_FILE_u, line=__LINE__) @@ -420,30 +474,38 @@ subroutine ice_import( importState, rc ) ! import ocn/ice fluxes - call state_getimport(importState, 'freezing_melting_potential', output=aflds, index=9, rc=rc) + call state_getimport(importState, 'freezing_melting_potential', output=aflds, index=9, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! import atm fluxes - call state_getimport(importState, 'mean_down_sw_vis_dir_flx', output=aflds, index=10, rc=rc) + call state_getimport(importState, 'mean_down_sw_vis_dir_flx', output=aflds, index=10, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_down_sw_ir_dir_flx', output=aflds, index=11, rc=rc) + call state_getimport(importState, 'mean_down_sw_ir_dir_flx', output=aflds, index=11, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_down_sw_vis_dif_flx', output=aflds, index=12, rc=rc) + call state_getimport(importState, 'mean_down_sw_vis_dif_flx', output=aflds, index=12, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_down_sw_ir_dif_flx', output=aflds, index=13, rc=rc) + call state_getimport(importState, 'mean_down_sw_ir_dif_flx', output=aflds, index=13, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_down_lw_flx', output=aflds, index=14, rc=rc) + call state_getimport(importState, 'mean_down_lw_flx', output=aflds, index=14, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_prec_rate', output=aflds, index=15, rc=rc) + call state_getimport(importState, 'mean_prec_rate', output=aflds, index=15, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'mean_fprec_rate', output=aflds, index=16, rc=rc) + call state_getimport(importState, 'mean_fprec_rate', output=aflds, index=16, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! perform a halo update @@ -479,7 +541,7 @@ subroutine ice_import( importState, rc ) end do !$OMP END PARALLEL DO - if ( State_fldChk(importState, 'Sa_ptem') .and. State_fldchk(importState,'air_density_height_lowest')) then + if ( State_fldChk(importState, 'Sa_ptem') .and. State_fldchk(importState,'air_density_height_lowest')) then !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1,ny_block @@ -509,7 +571,7 @@ subroutine ice_import( importState, rc ) endif end do !i end do !j - end do !iblk + end do !iblk !$OMP END PARALLEL DO end if @@ -568,34 +630,45 @@ subroutine ice_import( importState, rc ) ! bcphodry ungridded_index=2 ! bcphiwet ungridded_index=3 - ! bcphodry - call state_getimport(importState, 'Faxa_bcph', output=faero_atm, index=1, ungridded_index=2, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! bcphidry + bcphiwet - call state_getimport(importState, 'Faxa_bcph', output=faero_atm, index=2, ungridded_index=1, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_bcph', output=faero_atm, index=2, do_sum=.true., ungridded_index=3, rc=rc) + call state_getfldptr(importState, 'Faxa_bcph', fldptr=dataPtr2d, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + faero_atm(i,j,1,iblk) = dataPtr2d(2,n) * med2mod_areacor(n) ! bcphodry + faero_atm(i,j,2,iblk) = (dataptr2d(1,n) + dataPtr2d(3,n)) * med2mod_areacor(n) ! bcphidry + bcphiwet + end do + end do + end do end if ! Sum over all dry and wet dust fluxes from ath atmosphere if (State_FldChk(importState, 'Faxa_dstwet') .and. State_FldChk(importState, 'Faxa_dstdry')) then - call state_getimport(importState, 'Faxa_dstwet', output=faero_atm, index=3, ungridded_index=1, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstdry', output=faero_atm, index=3, do_sum=.true., ungridded_index=1, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstwet', output=faero_atm, index=3, do_sum=.true., ungridded_index=2, rc=rc) + call state_getfldptr(importState, 'Faxa_dstwet', fldptr=dataPtr2d_dstwet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstdry', output=faero_atm, index=3, do_sum=.true., ungridded_index=2, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstwet', output=faero_atm, index=3, do_sum=.true., ungridded_index=3, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstdry', output=faero_atm, index=3, do_sum=.true., ungridded_index=3, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstwet', output=faero_atm, index=3, do_sum=.true., ungridded_index=4, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Faxa_dstdry', output=faero_atm, index=3, do_sum=.true., ungridded_index=4, rc=rc) + call state_getfldptr(importState, 'Faxa_dstdry', fldptr=dataPtr2d_dstdry, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + faero_atm(i,j,3,iblk) = dataPtr2d_dstwet(1,n) + dataptr2d_dstdry(1,n) + & + dataPtr2d_dstwet(2,n) + dataptr2d_dstdry(2,n) + & + dataPtr2d_dstwet(3,n) + dataptr2d_dstdry(3,n) + & + dataPtr2d_dstwet(4,n) + dataptr2d_dstdry(4,n) + faero_atm(i,j,3,iblk) = faero_atm(i,j,3,iblk) * med2mod_areacor(n) + end do + end do + end do end if !------------------------------------------------------- @@ -614,12 +687,15 @@ subroutine ice_import( importState, rc ) call state_getimport(importState, 'inst_spec_humid_height_lowest_wiso', output=Qa_iso, index=3, ungridded_index=2, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return -! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=1, ungridded_index=3, rc=rc) -! if (ChkErr(rc,__LINE__,u_FILE_u)) return -! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=2, ungridded_index=1, rc=rc) -! if (ChkErr(rc,__LINE__,u_FILE_u)) return -! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=3, ungridded_index=2, rc=rc) -! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=1, ungridded_index=3, & + ! areacor=med2mod_areacor, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=2, ungridded_index=1, & + ! areacor=med2mod_areacor, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! call state_getimport(importState, 'mean_prec_rate_wiso', output=fiso_rain, index=3, ungridded_index=2, & + ! areacor=med2mod_areacor, rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_getimport(importState, 'mean_fprec_rate_wiso', output=fiso_atm, index=1, ungridded_index=3, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -628,11 +704,14 @@ subroutine ice_import( importState, rc ) call state_getimport(importState, 'mean_fprec_rate_wiso', output=fiso_atm, index=3, ungridded_index=2, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'So_roce_wiso', output=HDO_ocn , ungridded_index=3, rc=rc) + call state_getimport(importState, 'So_roce_wiso', output=HDO_ocn , ungridded_index=3, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'So_roce_wiso', output=H2_16O_ocn, ungridded_index=1, rc=rc) + call state_getimport(importState, 'So_roce_wiso', output=H2_16O_ocn, ungridded_index=1, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'So_roce_wiso', output=H2_18O_ocn, ungridded_index=2, rc=rc) + call state_getimport(importState, 'So_roce_wiso', output=H2_18O_ocn, ungridded_index=2, & + areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -682,7 +761,7 @@ subroutine ice_import( importState, rc ) #ifdef CESMCOUPLED ! Use shr_frz_mod for this Tf(:,:,iblk) = shr_frz_freezetemp(sss(:,:,iblk)) -#else +#else !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky) do iblk = 1, nblocks do j = 1,ny_block @@ -738,7 +817,6 @@ subroutine ice_import( importState, rc ) end subroutine ice_import !=============================================================================== - subroutine ice_export( exportState, rc ) ! input/output variables @@ -770,12 +848,13 @@ subroutine ice_export( exportState, rc ) if (io_dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) call icepack_query_parameters(Tffresh_out=Tffresh) -! call icepack_query_parameters(tfrz_option_out=tfrz_option, & -! modal_aero_out=modal_aero, z_tracers_out=z_tracers, skl_bgc_out=skl_bgc, & -! Tffresh_out=Tffresh) -! call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_iage_out=tr_iage, & -! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & -! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + ! call icepack_query_parameters(tfrz_option_out=tfrz_option, & + ! modal_aero_out=modal_aero, z_tracers_out=z_tracers, skl_bgc_out=skl_bgc, & + ! Tffresh_out=Tffresh) + ! call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_iage_out=tr_iage, & + ! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & + ! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=u_FILE_u, line=__LINE__) @@ -958,31 +1037,38 @@ subroutine ice_export( exportState, rc ) ! ------ ! Zonal air/ice stress - call state_setexport(exportState, 'stress_on_air_ice_zonal' , input=tauxa, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'stress_on_air_ice_zonal' , input=tauxa, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Meridional air/ice stress - call state_setexport(exportState, 'stress_on_air_ice_merid' , input=tauya, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'stress_on_air_ice_merid' , input=tauya, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Latent heat flux (atm into ice) - call state_setexport(exportState, 'mean_laten_heat_flx_atm_into_ice' , input=flat, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_laten_heat_flx_atm_into_ice' , input=flat, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Sensible heat flux (atm into ice) - call state_setexport(exportState, 'mean_sensi_heat_flx_atm_into_ice' , input=fsens, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sensi_heat_flx_atm_into_ice' , input=fsens, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! longwave outgoing (upward), average over ice fraction only - call state_setexport(exportState, 'mean_up_lw_flx_ice' , input=flwout, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_up_lw_flx_ice' , input=flwout, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Evaporative water flux (kg/m^2/s) - call state_setexport(exportState, 'mean_evap_rate_atm_into_ice' , input=evap, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_evap_rate_atm_into_ice' , input=evap, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Shortwave flux absorbed in ice and ocean (W/m^2) - call state_setexport(exportState, 'Faii_swnet' , input=fswabs, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'Faii_swnet' , input=fswabs, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! ------ @@ -990,43 +1076,53 @@ subroutine ice_export( exportState, rc ) ! ------ ! flux of shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn' , input=fswthru, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn' , input=fswthru, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! flux of vis dir shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn_vis_dir_flx' , input=fswthru_vdr, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn_vis_dir_flx' , input=fswthru_vdr, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! flux of vis dif shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn_vis_dif_flx' , input=fswthru_vdf, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn_vis_dif_flx' , input=fswthru_vdf, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! flux of ir dir shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn_ir_dir_flx' , input=fswthru_idr, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn_ir_dir_flx' , input=fswthru_idr, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! flux of ir dif shortwave through ice to ocean - call state_setexport(exportState, 'mean_sw_pen_to_ocn_ir_dif_flx' , input=fswthru_idf, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'mean_sw_pen_to_ocn_ir_dif_flx' , input=fswthru_idf, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! heat exchange with ocean - call state_setexport(exportState, 'net_heat_flx_to_ocn' , input=fhocn, lmask=tmask, ifrac=ailohi, rc=rc) + ! flux of heat exchange with ocean + call state_setexport(exportState, 'net_heat_flx_to_ocn' , input=fhocn, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! fresh water to ocean (h2o flux from melting) - call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate' , input=fresh, lmask=tmask, ifrac=ailohi, rc=rc) + ! flux fresh water to ocean (h2o flux from melting) + call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate' , input=fresh, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! salt to ocean (salt flux from melting) - call state_setexport(exportState, 'mean_salt_rate' , input=fsalt, lmask=tmask, ifrac=ailohi, rc=rc) + ! flux of salt to ocean (salt flux from melting) + call state_setexport(exportState, 'mean_salt_rate' , input=fsalt, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! stress n i/o zonal - call state_setexport(exportState, 'stress_on_ocn_ice_zonal' , input=tauxo, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'stress_on_ocn_ice_zonal' , input=tauxo, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! stress n i/o meridional - call state_setexport(exportState, 'stress_on_ocn_ice_merid' , input=tauyo, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'stress_on_ocn_ice_merid' , input=tauyo, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! ------ @@ -1035,19 +1131,22 @@ subroutine ice_export( exportState, rc ) ! hydrophobic bc if (State_FldChk(exportState, 'Fioi_bcpho')) then - call state_setexport(exportState, 'Fioi_bcpho' , input=faero_ocn, index=1, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'Fioi_bcpho' , input=faero_ocn, index=1, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if ! hydrophilic bc if (State_FldChk(exportState, 'Fioi_bcphi')) then - call state_setexport(exportState, 'Fioi_bcphi' , input=faero_ocn, index=2, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'Fioi_bcphi' , input=faero_ocn, index=2, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if ! dust if (State_FldChk(exportState, 'Fioi_flxdst')) then - call state_setexport(exportState, 'Fioi_flxdst' , input=faero_ocn, index=3, lmask=tmask, ifrac=ailohi, rc=rc) + call state_setexport(exportState, 'Fioi_flxdst' , input=faero_ocn, index=3, lmask=tmask, ifrac=ailohi, & + areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -1061,13 +1160,13 @@ subroutine ice_export( exportState, rc ) ! HDO => ungridded_index=3 call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate_wiso' , input=fiso_ocn, index=1, & - lmask=tmask, ifrac=ailohi, ungridded_index=3, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=3, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate_wiso' , input=fiso_ocn, index=2, & - lmask=tmask, ifrac=ailohi, ungridded_index=1, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=1, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_setexport(exportState, 'mean_fresh_water_to_ocean_rate_wiso' , input=fiso_ocn, index=3, & - lmask=tmask, ifrac=ailohi, ungridded_index=2, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=2, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -1078,16 +1177,16 @@ subroutine ice_export( exportState, rc ) if (State_FldChk(exportState, 'mean_evap_rate_atm_into_ice_wiso')) then ! Isotope evap to atm call state_setexport(exportState, 'mean_evap_rate_atm_into_ice_wiso' , input=fiso_evap, index=1, & - lmask=tmask, ifrac=ailohi, ungridded_index=3, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=3, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_setexport(exportState, 'mean_evap_rate_atm_into_ice_wiso' , input=fiso_evap, index=2, & - lmask=tmask, ifrac=ailohi, ungridded_index=1, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=1, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_setexport(exportState, 'mean_evap_rate_atm_into_ice_wiso' , input=fiso_evap, index=3, & - lmask=tmask, ifrac=ailohi, ungridded_index=2, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=2, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Isotope evap to atm + ! qref to atm call state_setexport(exportState, 'Si_qref_wiso' , input=Qref_iso, index=1, & lmask=tmask, ifrac=ailohi, ungridded_index=3, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -1115,7 +1214,7 @@ subroutine ice_export( exportState, rc ) ! Note: no need zero out pass-through fields over land for benefit of x2oacc fields in cpl hist files since ! the export state has been zeroed out at the beginning call state_setexport(exportState, 'mean_sw_pen_to_ocn_ifrac_n', input=fswthrun_ai, index=n, & - lmask=tmask, ifrac=ailohi, ungridded_index=n, rc=rc) + lmask=tmask, ifrac=ailohi, ungridded_index=n, areacor=mod2med_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end do end if @@ -1123,7 +1222,6 @@ subroutine ice_export( exportState, rc ) end subroutine ice_export !=============================================================================== - subroutine fldlist_add(num, fldlist, stdname, ungridded_lbound, ungridded_ubound) ! input/output variables @@ -1153,7 +1251,6 @@ subroutine fldlist_add(num, fldlist, stdname, ungridded_lbound, ungridded_ubound end subroutine fldlist_add !=============================================================================== - subroutine fldlist_realize(state, fldList, numflds, flds_scalar_name, flds_scalar_num, mesh, grid, tag, rc) use NUOPC, only : NUOPC_IsConnected, NUOPC_Realize @@ -1284,7 +1381,6 @@ end subroutine SetScalarField end subroutine fldlist_realize !=============================================================================== - logical function State_FldChk(State, fldname) ! ---------------------------------------------- ! Determine if field is in state @@ -1299,27 +1395,25 @@ logical function State_FldChk(State, fldname) ! ---------------------------------------------- call ESMF_StateGet(State, trim(fldname), itemType) - State_FldChk = (itemType /= ESMF_STATEITEM_NOTFOUND) end function State_FldChk !=============================================================================== - - subroutine state_getimport_4d_output(state, fldname, output, index, do_sum, ungridded_index, rc) + subroutine state_getimport_4d_output(state, fldname, output, index, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map import state field to output array ! ---------------------------------------------- ! input/output variables - type(ESMF_State) , intent(in) :: state - character(len=*) , intent(in) :: fldname - real (kind=dbl_kind) , intent(inout) :: output(:,:,:,:) - integer , intent(in) :: index - logical, optional , intent(in) :: do_sum - integer, optional , intent(in) :: ungridded_index - integer , intent(out) :: rc + type(ESMF_State) , intent(in) :: state + character(len=*) , intent(in) :: fldname + real (kind=dbl_kind) , intent(inout) :: output(:,:,:,:) + integer , intent(in) :: index + integer, optional , intent(in) :: ungridded_index + real(kind=dbl_kind), optional , intent(in) :: areacor(:) + integer , intent(out) :: rc ! local variables type(block) :: this_block ! block information for current block @@ -1327,9 +1421,7 @@ subroutine state_getimport_4d_output(state, fldname, output, index, do_sum, ungr integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - real(kind=dbl_kind), pointer :: dataPtr3d(:,:,:) ! grid - real(kind=dbl_kind), pointer :: dataPtr4d(:,:,:,:) ! grid - character(len=*), parameter :: subname='(ice_import_export:state_getimport)' + character(len=*), parameter :: subname='(ice_import_export:state_getimport_4d_output)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1337,103 +1429,62 @@ subroutine state_getimport_4d_output(state, fldname, output, index, do_sum, ungr ! check that fieldname exists if (.not. State_FldChk(state, trim(fldname))) return - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataPtr2d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataPtr1d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + ! get field pointer + if (present(ungridded_index)) then + call state_getfldptr(state, trim(fldname), dataPtr2d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call state_getfldptr(state, trim(fldname), dataPtr1d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - ! set values of output array - n=0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - if (present(do_sum)) then ! do sum - if (present(ungridded_index)) then - output(i,j,index,iblk) = output(i,j,index, iblk) + dataPtr2d(ungridded_index,n) - else - output(i,j,index,iblk) = output(i,j,index, iblk) + dataPtr1d(n) - end if - else ! do not do sum - if (present(ungridded_index)) then - output(i,j,index,iblk) = dataPtr2d(ungridded_index,n) - else - output(i,j,index,iblk) = dataPtr1d(n) - end if - end if - end do + ! set values of output array + n=0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + if (present(ungridded_index)) then + output(i,j,index,iblk) = dataPtr2d(ungridded_index,n) + else + output(i,j,index,iblk) = dataPtr1d(n) + end if end do end do - - else if (geomtype == ESMF_GEOMTYPE_GRID) then - - ! get field pointer + end do + if (present(areacor)) then if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataptr4d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + do n = 1,size(dataptr2d,dim=2) + dataPtr2d(:,n) = dataPtr2d(:,n) * areacor(n) + end do else - call state_getfldptr(state, trim(fldname), dataptr3d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - ! set values of output array - do iblk = 1,nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo,jhi - do i = ilo,ihi - i1 = i - ilo + 1 - j1 = j - jlo + 1 - if (present(do_sum)) then - if (present(ungridded_index)) then - output(i,j,index,iblk) = output(i,j,index,iblk) + dataPtr4d(i1,j1,iblk,ungridded_index) - else - output(i,j,index,iblk) = output(i,j,index,iblk) + dataPtr3d(i1,j1,iblk) - end if - else - if (present(ungridded_index)) then - output(i,j,index,iblk) = dataPtr4d(i1,j1,iblk,ungridded_index) - else - output(i,j,index,iblk) = dataPtr3d(i1,j1,iblk) - end if - end if - end do + do n = 1,size(dataptr1d) + dataPtr1d(n) = dataPtr1d(n) * areacor(n) end do - end do - + end if end if end subroutine state_getimport_4d_output !=============================================================================== - - subroutine state_getimport_3d_output(state, fldname, output, do_sum, ungridded_index, rc) + subroutine state_getimport_3d_output(state, fldname, output, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map import state field to output array ! ---------------------------------------------- ! input/output variables - type(ESMF_State) , intent(in) :: state - character(len=*) , intent(in) :: fldname - real (kind=dbl_kind) , intent(inout) :: output(:,:,:) - logical, optional , intent(in) :: do_sum - integer, optional , intent(in) :: ungridded_index - integer , intent(out) :: rc + type(ESMF_State) , intent(in) :: state + character(len=*) , intent(in) :: fldname + real (kind=dbl_kind) , intent(inout) :: output(:,:,:) + integer, optional , intent(in) :: ungridded_index + real(kind=dbl_kind), optional , intent(in) :: areacor(:) + integer , intent(out) :: rc ! local variables type(block) :: this_block ! block information for current block @@ -1441,8 +1492,6 @@ subroutine state_getimport_3d_output(state, fldname, output, do_sum, ungridded_i integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - real(kind=dbl_kind), pointer :: dataPtr3d(:,:,:) ! grid - real(kind=dbl_kind), pointer :: dataPtr4d(:,:,:,:) ! grid character(len=*) , parameter :: subname='(ice_import_export:state_getimport)' ! ---------------------------------------------- @@ -1451,83 +1500,53 @@ subroutine state_getimport_3d_output(state, fldname, output, do_sum, ungridded_i ! check that fieldname exists if (.not. State_FldChk(state, trim(fldname))) return - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataPtr2d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataPtr1d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + ! get field pointer + if (present(ungridded_index)) then + call state_getfldptr(state, trim(fldname), dataPtr2d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call state_getfldptr(state, trim(fldname), dataPtr1d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - ! determine output array - n=0 + ! determine output array + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + if (present(ungridded_index)) then + output(i,j,iblk) = dataPtr2d(ungridded_index,n) + else + output(i,j,iblk) = dataPtr1d(n) + end if + end do + end do + end do + if (present(areacor)) then + n = 0 do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi do j = jlo, jhi do i = ilo, ihi - n = n+1 - if (present(do_sum) .and. present(ungridded_index)) then - output(i,j,iblk) = output(i,j,iblk) + dataPtr2d(ungridded_index,n) - else if (present(do_sum)) then - output(i,j,iblk) = output(i,j,iblk) + dataPtr1d(n) - else if (present(ungridded_index)) then - output(i,j,iblk) = dataPtr2d(ungridded_index,n) - else - output(i,j,iblk) = dataPtr1d(n) - end if + n = n + 1 + output(i,j,iblk) = output(i,j,iblk) * areacor(n) end do end do end do - - else if (geomtype == ESMF_GEOMTYPE_GRID) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataptr4d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataptr3d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - ! set values of output array - do iblk = 1,nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo,jhi - do i = ilo,ihi - i1 = i - ilo + 1 - j1 = j - jlo + 1 - if (present(do_sum) .and. present(ungridded_index)) then - output(i,j,iblk) = output(i,j,iblk) + dataPtr4d(i1,j1,iblk,ungridded_index) - else if (present(do_sum)) then - output(i,j,iblk) = output(i,j,iblk) + dataPtr3d(i1,j1,iblk) - else if (present(ungridded_index)) then - output(i,j,iblk) = dataPtr4d(i1,j1,iblk, ungridded_index) - else - output(i,j,iblk) = dataPtr3d(i1,j1,iblk) - end if - end do - end do - end do - end if end subroutine state_getimport_3d_output !=============================================================================== - - subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, ungridded_index, rc) + subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map 4d input array to export state field @@ -1541,6 +1560,7 @@ subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, logical , optional, intent(in) :: lmask(:,:,:) real(kind=dbl_kind) , optional, intent(in) :: ifrac(:,:,:) integer , optional, intent(in) :: ungridded_index + real(kind=dbl_kind) , optional, intent(in) :: areacor(:) integer , intent(out) :: rc ! local variables @@ -1549,8 +1569,6 @@ subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, integer :: i, j, iblk, n, i1, j1 ! indices real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - real(kind=dbl_kind), pointer :: dataPtr3d(:,:,:) ! grid - real(kind=dbl_kind), pointer :: dataPtr4d(:,:,:,:) ! grid character(len=*), parameter :: subname='(ice_import_export:state_setexport)' ! ---------------------------------------------- @@ -1559,93 +1577,60 @@ subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, ! check that fieldname exists if (.not. State_FldChk(state, trim(fldname))) return - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataPtr2d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataPtr1d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + ! get field pointer + if (present(ungridded_index)) then + call state_getfldptr(state, trim(fldname), dataPtr2d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call state_getfldptr(state, trim(fldname), dataPtr1d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - ! set values of field pointer - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - if (present(lmask) .and. present(ifrac)) then - if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then - if (present(ungridded_index)) then - dataPtr2d(ungridded_index,n) = input(i,j,index,iblk) - else - dataPtr1d(n) = input(i,j,index,iblk) - end if - end if - else + ! set values of field pointer + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + if (present(lmask) .and. present(ifrac)) then + if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then if (present(ungridded_index)) then dataPtr2d(ungridded_index,n) = input(i,j,index,iblk) else dataPtr1d(n) = input(i,j,index,iblk) end if end if - end do + else + if (present(ungridded_index)) then + dataPtr2d(ungridded_index,n) = input(i,j,index,iblk) + else + dataPtr1d(n) = input(i,j,index,iblk) + end if + end if end do end do - - else if (geomtype == ESMF_GEOMTYPE_GRID) then - + end do + if (present(areacor)) then if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataptr4d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + do n = 1,size(dataptr2d,dim=2) + dataPtr2d(:,n) = dataPtr2d(:,n) * areacor(n) + end do else - call state_getfldptr(state, trim(fldname), dataptr3d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - do iblk = 1,nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo,jhi - do i = ilo,ihi - i1 = i - ilo + 1 - j1 = j - jlo + 1 - if (present(lmask) .and. present(ifrac)) then - if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then - if (present(ungridded_index)) then - dataPtr4d(i1,j1,iblk,ungridded_index) = input(i,j,index,iblk) - end if - else - dataPtr3d(i1,j1,iblk) = input(i,j,index,iblk) - end if - else - if (present(ungridded_index)) then - dataPtr4d(i1,j1,iblk,ungridded_index) = input(i,j,index,iblk) - else - dataPtr3d(i1,j1,iblk) = input(i,j,index,iblk) - end if - end if - end do + do n = 1,size(dataptr1d) + dataPtr1d(n) = dataPtr1d(n) * areacor(n) end do - end do - + end if end if end subroutine state_setexport_4d_input !=============================================================================== - - subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridded_index, rc) + subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map 3d input array to export state field @@ -1658,6 +1643,7 @@ subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridd logical , optional , intent(in) :: lmask(:,:,:) real(kind=dbl_kind) , optional , intent(in) :: ifrac(:,:,:) integer , optional , intent(in) :: ungridded_index + real(kind=dbl_kind) , optional , intent(in) :: areacor(:) integer , intent(out) :: rc ! local variables @@ -1666,8 +1652,6 @@ subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridd integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - real(kind=dbl_kind), pointer :: dataPtr3d(:,:,:) ! grid - real(kind=dbl_kind), pointer :: dataPtr4d(:,:,:,:) ! grid character(len=*), parameter :: subname='(ice_import_export:state_setexport)' ! ---------------------------------------------- @@ -1676,92 +1660,58 @@ subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridd ! check that fieldname exists if (.not. State_FldChk(state, trim(fldname))) return - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! get field pointer - if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataPtr2d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - else - call state_getfldptr(state, trim(fldname), dataPtr1d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + ! get field pointer + if (present(ungridded_index)) then + call state_getfldptr(state, trim(fldname), dataPtr2d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call state_getfldptr(state, trim(fldname), dataPtr1d, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - if (present(lmask) .and. present(ifrac)) then - if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then - if (present(ungridded_index)) then - dataPtr2d(ungridded_index,n) = input(i,j,iblk) - else - dataPtr1d(n) = input(i,j,iblk) - end if - end if - else + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + if (present(lmask) .and. present(ifrac)) then + if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then if (present(ungridded_index)) then dataPtr2d(ungridded_index,n) = input(i,j,iblk) else dataPtr1d(n) = input(i,j,iblk) end if end if - end do + else + if (present(ungridded_index)) then + dataPtr2d(ungridded_index,n) = input(i,j,iblk) + else + dataPtr1d(n) = input(i,j,iblk) + end if + end if end do end do - - else if (geomtype == ESMF_GEOMTYPE_GRID) then - - ! get field pointer + end do + if (present(areacor)) then if (present(ungridded_index)) then - call state_getfldptr(state, trim(fldname), dataptr4d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + do n = 1,size(dataptr2d,dim=2) + dataPtr2d(:,n) = dataPtr2d(:,n) * areacor(n) + end do else - call state_getfldptr(state, trim(fldname), dataptr3d, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - do iblk = 1,nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - do j = jlo,jhi - do i = ilo,ihi - i1 = i - ilo + 1 - j1 = j - jlo + 1 - if (present(lmask) .and. present(ifrac)) then - if ( lmask(i,j,iblk) .and. ifrac(i,j,iblk) > c0 ) then - if (present(ungridded_index)) then - dataPtr4d(i1,j1,iblk,ungridded_index) = input(i,j,iblk) - else - dataPtr3d(i1,j1,iblk) = input(i,j,iblk) - end if - end if - else - if (present(ungridded_index)) then - dataPtr4d(i1,j1,iblk,ungridded_index) = input(i,j,iblk) - else - dataPtr3d(i1,j1,iblk) = input(i,j,iblk) - end if - end if - end do + do n = 1,size(dataptr1d) + dataPtr1d(n) = dataPtr1d(n) * areacor(n) end do - end do - + end if end if end subroutine state_setexport_3d_input !=============================================================================== - subroutine State_GetFldPtr_1d(State, fldname, fldptr, rc) ! ---------------------------------------------- ! Get pointer to a state field @@ -1785,10 +1735,10 @@ subroutine State_GetFldPtr_1d(State, fldname, fldptr, rc) call ESMF_FieldGet(lfield, farrayPtr=fldptr, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + end subroutine State_GetFldPtr_1d !=============================================================================== - subroutine State_GetFldPtr_2d(State, fldname, fldptr, rc) ! ---------------------------------------------- ! Get pointer to a state field @@ -1807,65 +1757,12 @@ subroutine State_GetFldPtr_2d(State, fldname, fldptr, rc) rc = ESMF_SUCCESS - call ESMF_StateGet(State, itemName=trim(fldname), field=lfield, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call ESMF_FieldGet(lfield, farrayPtr=fldptr, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end subroutine State_GetFldPtr_2d - - !=============================================================================== - - subroutine State_GetFldPtr_3d(State, fldname, fldptr, rc) - ! ---------------------------------------------- - ! Get pointer to a state field - ! ---------------------------------------------- - - ! input/output variables - type(ESMF_State) , intent(in) :: State - character(len=*) , intent(in) :: fldname - real(kind=dbl_kind) , pointer , intent(inout) :: fldptr(:,:,:) - integer , optional , intent(out) :: rc - - ! local variables - type(ESMF_Field) :: lfield - character(len=*),parameter :: subname='(ice_import_export:State_GetFldPtr_3d)' - ! ---------------------------------------------- - - rc = ESMF_SUCCESS - call ESMF_StateGet(State, itemName=trim(fldname), field=lfield, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldGet(lfield, farrayPtr=fldptr, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - end subroutine State_GetFldPtr_3d - !=============================================================================== - - subroutine State_GetFldPtr_4d(State, fldname, fldptr, rc) - ! ---------------------------------------------- - ! Get pointer to a state field - ! ---------------------------------------------- - - ! input/output variables - type(ESMF_State) , intent(in) :: State - character(len=*) , intent(in) :: fldname - real(kind=dbl_kind) , pointer , intent(inout) :: fldptr(:,:,:,:) - integer , optional , intent(out) :: rc - - ! local variables - type(ESMF_Field) :: lfield - character(len=*),parameter :: subname='(ice_import_export:State_GetFldPtr_3d)' - ! ---------------------------------------------- - - rc = ESMF_SUCCESS - - call ESMF_StateGet(State, itemName=trim(fldname), field=lfield, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call ESMF_FieldGet(lfield, farrayPtr=fldptr, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end subroutine State_GetFldPtr_4d + end subroutine State_GetFldPtr_2d end module ice_import_export diff --git a/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 b/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 index b02865320..b9dfda5c7 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 @@ -1,14 +1,17 @@ module ice_mesh_mod use ESMF - use ice_kinds_mod , only : dbl_kind, int_kind + use NUOPC , only : NUOPC_CompAttributeGet + use ice_kinds_mod , only : dbl_kind, int_kind, char_len, char_len_long use ice_domain_size , only : nx_global, ny_global, max_blocks use ice_domain , only : nblocks, blocks_ice, distrb_info use ice_blocks , only : block, get_block, nx_block, ny_block, nblocks_x, nblocks_y use ice_shr_methods , only : chkerr use ice_fileunits , only : nu_diag use ice_communicate , only : my_task, master_task - use ice_initfc , only : icepack_query_parameters + use ice_exit , only : abort_ice + use icepack_intfc , only : icepack_query_parameters + use icepack_intfc , only : icepack_warnings_flush, icepack_warnings_aborted implicit none private @@ -200,7 +203,7 @@ subroutine ice_mesh_set_distgrid(localpet, npes, distgrid, rc) end subroutine ice_mesh_set_distgrid !======================================================================= - subroutine ice_mesh_setmask_from_maskfile(ice_mesh, ice_maskfile, rc) + subroutine ice_mesh_setmask_from_maskfile(gcomp, ice_mesh, rc) use ice_scam , only : scmlat, scmlon, single_column use ice_grid , only : tlon, tlat, hm, tarea, ULON, ULAT, HTN, HTE, ANGLE, ANGLET @@ -209,18 +212,18 @@ subroutine ice_mesh_setmask_from_maskfile(ice_mesh, ice_maskfile, rc) use ice_grid , only : kmt_file, makemask, tmask use ice_boundary , only : ice_HaloUpdate use ice_domain , only : blocks_ice, nblocks, halo_info, distrb_info - use ice_constants , only : c0, c1, c2, p5, radius + use ice_constants , only : c0, c1, c2, p25, radius use ice_constants , only : field_loc_center, field_type_scalar use ice_read_write , only : ice_open_nc, ice_close_nc - use ice_exit , only : abort_ice use netcdf ! input/output variables - type(ESMF_Mesh) , intent(in) :: ice_mesh - character(len=*) , intent(in) :: ice_maskfile - integer , intent(out) :: rc + type(ESMF_GridComp) , intent(inout) :: gcomp + type(ESMF_Mesh) , intent(in) :: ice_mesh + integer , intent(out) :: rc ! local variables + character(len=char_len_long) :: ice_maskfile integer :: i,j,n integer :: iblk, jblk ! indices integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain @@ -250,13 +253,41 @@ subroutine ice_mesh_setmask_from_maskfile(ice_mesh, ice_maskfile, rc) integer (int_kind), pointer :: ice_mask(:) real(dbl_kind) , pointer :: ice_frac(:) real(dbl_kind) :: pi + real(dbl_kind) :: c180 real(dbl_kind) :: puny real(dbl_kind) :: deg_to_rad + logical :: isPresent, isSet + character(len=char_len_long) :: cvalue character(len=*), parameter :: subname = ' ice_mesh_setmask_from_maskfile' !--------------------------------------------------- rc = ESMF_SUCCESS + ! Determine mask input file and create the mask + call NUOPC_CompAttributeGet(gcomp, name='mesh_mask', value=ice_maskfile, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (my_task == master_task) then + write(nu_diag,*)'mask file for cice domain is ',trim(ice_maskfile) + end if + + ! Determine if single column + call NUOPC_CompAttributeGet(gcomp, name='single_column', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) single_column + if (single_column) then + call NUOPC_CompAttributeGet(gcomp, name='scmlon', value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) scmlon + call NUOPC_CompAttributeGet(gcomp, name='scmlat', value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) scmlat + end if + else + single_column = .false. + end if + ! Determine start/count to read in for either single column or global lat-lon grid ! If single_column, then assume that only master_task is used since there is only one task @@ -288,9 +319,10 @@ subroutine ice_mesh_setmask_from_maskfile(ice_mesh, ice_maskfile, rc) allocate(ocn_gridcell_frac(nx_block,ny_block,max_blocks)) ! Get required constants - call icepack_query_parameters(pi_out=pi, pi2_out=pi2, puny_out=puny, deg_to_rad_out=deg_to_rad) + call icepack_query_parameters(pi_out=pi, puny_out=puny, c180_out=c180) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) + deg_to_rad = pi/c180 ! Set tlon, tlat, tarea, hm and ocn_gridcell_frac ! Convert mesh areas from radians^2 to m^2 (tarea is in m^2) @@ -641,6 +673,7 @@ subroutine ice_mesh_check(ice_mesh, rc) real(dbl_kind) :: diff_lon real(dbl_kind) :: diff_lat real(dbl_kind) :: rad_to_deg + character(len=*), parameter :: subname = ' ice_mesh_check: ' !--------------------------------------------------- ! error check differences between internally generated lons and those read in @@ -680,11 +713,13 @@ subroutine ice_mesh_check(ice_mesh, rc) if ( (diff_lon > 1.e2 .and. abs(diff_lon - 360.) > 1.e-1) .or.& (diff_lon > 1.e-3 .and. diff_lon < c1) ) then write(6,100)n,lonMesh(n),lon(n), diff_lon - !call shr_sys_abort() + call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) end if if (abs(latMesh(n) - lat(n)) > 1.e-1) then write(6,101)n,latMesh(n),lat(n), abs(latMesh(n)-lat(n)) - !call shr_sys_abort() + call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) end if enddo From 0f90ce174461fe9e3e64474faf9ef9465292b0cf Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Mon, 1 Feb 2021 09:43:37 -0700 Subject: [PATCH 4/4] more cleanup --- .../drivers/nuopc/cmeps/ice_import_export.F90 | 67 +++++++++++-------- 1 file changed, 38 insertions(+), 29 deletions(-) diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index 6e584eb90..82afad28d 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -59,14 +59,14 @@ module ice_import_export private :: state_getfldptr interface state_getimport - module procedure state_getimport_4d_output - module procedure state_getimport_3d_output + module procedure state_getimport_4d + module procedure state_getimport_3d end interface state_getimport private :: state_getimport interface state_setexport - module procedure state_setexport_4d_input - module procedure state_setexport_3d_input + module procedure state_setexport_4d + module procedure state_setexport_3d end interface state_setexport private :: state_setexport @@ -291,6 +291,7 @@ subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc real(dbl_kind), allocatable :: mesh_areas(:) real(dbl_kind), allocatable :: model_areas(:) real(dbl_kind), pointer :: dataptr(:) + integer :: num_ice character(len=*), parameter :: subname='(ice_import_export:realize_fields)' !--------------------------------------------------------------------------- @@ -355,9 +356,10 @@ subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc enddo enddo enddo + num_ice = n ! Determine flux correction factors (module variables) - do n = 1,numOwnedElements + do n = 1,num_ice if (model_areas(n) == mesh_areas(n)) then mod2med_areacor(n) = c1 med2mod_areacor(n) = c1 @@ -1400,7 +1402,7 @@ logical function State_FldChk(State, fldname) end function State_FldChk !=============================================================================== - subroutine state_getimport_4d_output(state, fldname, output, index, ungridded_index, areacor, rc) + subroutine state_getimport_4d(state, fldname, output, index, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map import state field to output array @@ -1421,7 +1423,7 @@ subroutine state_getimport_4d_output(state, fldname, output, index, ungridded_in integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - character(len=*), parameter :: subname='(ice_import_export:state_getimport_4d_output)' + character(len=*), parameter :: subname='(ice_import_export:state_getimport_4d)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1458,21 +1460,24 @@ subroutine state_getimport_4d_output(state, fldname, output, index, ungridded_in end do end do if (present(areacor)) then - if (present(ungridded_index)) then - do n = 1,size(dataptr2d,dim=2) - dataPtr2d(:,n) = dataPtr2d(:,n) * areacor(n) - end do - else - do n = 1,size(dataptr1d) - dataPtr1d(n) = dataPtr1d(n) * areacor(n) + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n + 1 + output(i,j,index,iblk) = output(i,j,index,iblk) * areacor(n) + end do end do - end if + end do end if - end subroutine state_getimport_4d_output + end subroutine state_getimport_4d !=============================================================================== - subroutine state_getimport_3d_output(state, fldname, output, ungridded_index, areacor, rc) + subroutine state_getimport_3d(state, fldname, output, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map import state field to output array @@ -1492,7 +1497,7 @@ subroutine state_getimport_3d_output(state, fldname, output, ungridded_index, ar integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - character(len=*) , parameter :: subname='(ice_import_export:state_getimport)' + character(len=*) , parameter :: subname='(ice_import_export:state_getimport_3d)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1543,10 +1548,10 @@ subroutine state_getimport_3d_output(state, fldname, output, ungridded_index, ar end do end if - end subroutine state_getimport_3d_output + end subroutine state_getimport_3d !=============================================================================== - subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, ungridded_index, areacor, rc) + subroutine state_setexport_4d(state, fldname, input, index, lmask, ifrac, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map 4d input array to export state field @@ -1569,7 +1574,8 @@ subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, integer :: i, j, iblk, n, i1, j1 ! indices real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - character(len=*), parameter :: subname='(ice_import_export:state_setexport)' + integer :: num_ice + character(len=*), parameter :: subname='(ice_import_export:state_setexport_4d)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1616,21 +1622,22 @@ subroutine state_setexport_4d_input(state, fldname, input, index, lmask, ifrac, end do end do if (present(areacor)) then + num_ice = n if (present(ungridded_index)) then - do n = 1,size(dataptr2d,dim=2) + do n = 1,num_ice dataPtr2d(:,n) = dataPtr2d(:,n) * areacor(n) end do else - do n = 1,size(dataptr1d) + do n = 1,num_ice dataPtr1d(n) = dataPtr1d(n) * areacor(n) end do end if end if - end subroutine state_setexport_4d_input + end subroutine state_setexport_4d !=============================================================================== - subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridded_index, areacor, rc) + subroutine state_setexport_3d(state, fldname, input, lmask, ifrac, ungridded_index, areacor, rc) ! ---------------------------------------------- ! Map 3d input array to export state field @@ -1652,7 +1659,8 @@ subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridd integer :: i, j, iblk, n, i1, j1 ! incides real(kind=dbl_kind), pointer :: dataPtr1d(:) ! mesh real(kind=dbl_kind), pointer :: dataPtr2d(:,:) ! mesh - character(len=*), parameter :: subname='(ice_import_export:state_setexport)' + integer :: num_ice + character(len=*), parameter :: subname='(ice_import_export:state_setexport_3d)' ! ---------------------------------------------- rc = ESMF_SUCCESS @@ -1697,19 +1705,20 @@ subroutine state_setexport_3d_input(state, fldname, input, lmask, ifrac, ungridd end do end do end do + num_ice = n if (present(areacor)) then if (present(ungridded_index)) then - do n = 1,size(dataptr2d,dim=2) + do n = 1,num_ice dataPtr2d(:,n) = dataPtr2d(:,n) * areacor(n) end do else - do n = 1,size(dataptr1d) + do n = 1,num_ice dataPtr1d(n) = dataPtr1d(n) * areacor(n) end do end if end if - end subroutine state_setexport_3d_input + end subroutine state_setexport_3d !=============================================================================== subroutine State_GetFldPtr_1d(State, fldname, fldptr, rc)