diff --git a/cicecore/cicedyn/dynamics/ice_dyn_evp1d.F90 b/cicecore/cicedyn/dynamics/ice_dyn_evp1d.F90 index 223ef2849..18316640e 100644 --- a/cicecore/cicedyn/dynamics/ice_dyn_evp1d.F90 +++ b/cicecore/cicedyn/dynamics/ice_dyn_evp1d.F90 @@ -94,7 +94,7 @@ subroutine dyn_evp1d_init endif ! gather from blks to global - call gather_static(G_uarear, G_dxT, G_dyT, G_tmask) + call gather_static(G_uarear, G_dxT, G_dyT, G_tmask) ! calculate number of water points (T and U). Only needed for the static version ! tmask in ocean/ice @@ -349,7 +349,6 @@ subroutine evp1d_alloc_static_na(na0) call abort_ice(subname//' ERROR: allocating', file=__FILE__, line=__LINE__) endif - allocate(indxTi(1:na0), & indxTj(1:na0), & stat=ierr) @@ -628,6 +627,11 @@ subroutine gather_static(G_uarear, G_dxT, G_dyT, G_Tmask) character(len=*), parameter :: subname = '(gather_static)' + G_uarear = c0 + G_dyT = c0 + G_dxT = c0 + G_tmask = .false. + ! copy from distributed I_* to G_* call gather_global_ext(G_uarear, uarear, master_task, distrb_info) call gather_global_ext(G_dxT , dxT , master_task, distrb_info) @@ -977,6 +981,37 @@ subroutine convert_1d_2d_dyn(na0 , navel0 , integer(kind=int_kind) :: lo, up, iw, i, j character(len=*), parameter :: subname = '(convert_1d_2d_dyn)' + G_stressp_1 = c0 + G_stressp_2 = c0 + G_stressp_3 = c0 + G_stressp_4 = c0 + G_stressm_1 = c0 + G_stressm_2 = c0 + G_stressm_3 = c0 + G_stressm_4 = c0 + G_stress12_1 = c0 + G_stress12_2 = c0 + G_stress12_3 = c0 + G_stress12_4 = c0 + G_strength = c0 + G_cdn_ocn = c0 + G_aiu = c0 + G_uocn = c0 + G_vocn = c0 + G_waterxU = c0 + G_wateryU = c0 + G_forcexU = c0 + G_forceyU = c0 + G_umassdti = c0 + G_fmU = c0 + G_strintxU = c0 + G_strintyU = c0 + G_Tbu = c0 + G_uvel = c0 + G_vvel = c0 + G_taubxU = c0 + G_taubyU = c0 + lo=1 up=na0 do iw = lo, up diff --git a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 index 40f49877d..81d603124 100644 --- a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 @@ -193,6 +193,13 @@ subroutine alloc_dyn_shared stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of memory') + uvel_init = c0 + vvel_init = c0 + iceTmask = .false. + iceUmask = .false. + fcor_blk = c0 + DminTarea = c0 + allocate( & fld2(nx_block,ny_block,2,max_blocks), & fld3(nx_block,ny_block,3,max_blocks), & @@ -200,6 +207,10 @@ subroutine alloc_dyn_shared stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of memory') + fld2 = c0 + fld3 = c0 + fld4 = c0 + allocate( & cyp(nx_block,ny_block,max_blocks), & ! 1.5*HTE - 0.5*HTW cxp(nx_block,ny_block,max_blocks), & ! 1.5*HTN - 0.5*HTS @@ -208,12 +219,19 @@ subroutine alloc_dyn_shared stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of memory') + cyp = c0 + cxp = c0 + cym = c0 + cxm = c0 + if (grid_ice == 'B' .and. evp_algorithm == "standard_2d") then allocate( & dxhy(nx_block,ny_block,max_blocks), & ! 0.5*(HTE - HTW) dyhx(nx_block,ny_block,max_blocks), & ! 0.5*(HTN - HTS) stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of memory') + dxhy = c0 + dyhx = c0 endif if (grid_ice == 'CD' .or. grid_ice == 'C') then @@ -228,6 +246,14 @@ subroutine alloc_dyn_shared fcorN_blk (nx_block,ny_block,max_blocks), & ! Coriolis stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of memory') + uvelE_init = c0 + vvelE_init = c0 + uvelN_init = c0 + vvelN_init = c0 + iceEmask = .false. + iceNmask = .false. + fcorE_blk = c0 + fcorN_blk = c0 endif end subroutine alloc_dyn_shared diff --git a/cicecore/cicedyn/general/ice_flux.F90 b/cicecore/cicedyn/general/ice_flux.F90 index ff71a4a4d..16d3a6edc 100644 --- a/cicecore/cicedyn/general/ice_flux.F90 +++ b/cicecore/cicedyn/general/ice_flux.F90 @@ -625,40 +625,270 @@ subroutine alloc_flux stat=ierr) if (ierr/=0) call abort_ice('(alloc_flux): Out of memory') - if (grid_ice == "CD" .or. grid_ice == "C") & + strax = c0 + stray = c0 + uocn = c0 + vocn = c0 + ss_tltx = c0 + ss_tlty = c0 + hwater = c0 + strairxT = c0 + strairyT = c0 + strocnxT_iavg= c0 + strocnyT_iavg= c0 + sig1 = c0 + sig2 = c0 + sigP = c0 + taubxU = c0 + taubyU = c0 + strairxU = c0 + strairyU = c0 + strocnxU = c0 + strocnyU = c0 + strtltxU = c0 + strtltyU = c0 + strintxU = c0 + strintyU = c0 + daidtd = c0 + dvidtd = c0 + dvsdtd = c0 + dagedtd = c0 + dardg1dt = c0 + dardg2dt = c0 + dvirdgdt = c0 + opening = c0 + stressp_1 = c0 + stressp_2 = c0 + stressp_3 = c0 + stressp_4 = c0 + stressm_1 = c0 + stressm_2 = c0 + stressm_3 = c0 + stressm_4 = c0 + stress12_1 = c0 + stress12_2 = c0 + stress12_3 = c0 + stress12_4 = c0 + fmU = c0 + TbU = c0 + zlvl = c0 + zlvs = c0 + uatm = c0 + vatm = c0 + wind = c0 + potT = c0 + Tair = c0 + Qa = c0 + rhoa = c0 + swvdr = c0 + swvdf = c0 + swidr = c0 + swidf = c0 + swuvrdr = c0 + swuvrdf = c0 + swpardr = c0 + swpardf = c0 + flw = c0 + frain = c0 + fsnow = c0 + sss = c0 + sst = c0 + frzmlt = c0 + frzmlt_init= c0 + Tf = c0 + qdp = c0 + hmix = c0 + daice_da = c0 + fsens = c0 + flat = c0 + fswabs = c0 + fswint_ai = c0 + flwout = c0 + Tref = c0 + Qref = c0 + Uref = c0 + evap = c0 + evaps = c0 + evapi = c0 + alvdr = c0 + alidr = c0 + alvdf = c0 + alidf = c0 + alvdr_ai = c0 + alidr_ai = c0 + alvdf_ai = c0 + alidf_ai = c0 + albice = c0 + albsno = c0 + albpnd = c0 + apeff_ai = c0 + snowfrac = c0 + alvdr_init = c0 + alidr_init = c0 + alvdf_init = c0 + alidf_init = c0 + fpond = c0 + fresh = c0 + fsalt = c0 + fhocn = c0 + fsloss = c0 + fswthru = c0 + fswthru_vdr= c0 + fswthru_vdf= c0 + fswthru_idr= c0 + fswthru_idf= c0 + fswthru_uvrdr = c0 + fswthru_uvrdf = c0 + fswthru_pardr = c0 + fswthru_pardf = c0 + scale_factor = c0 + strairx_ocn= c0 + strairy_ocn= c0 + fsens_ocn = c0 + flat_ocn = c0 + flwout_ocn = c0 + evap_ocn = c0 + alvdr_ocn = c0 + alidr_ocn = c0 + alvdf_ocn = c0 + alidf_ocn = c0 + Tref_ocn = c0 + Qref_ocn = c0 + fsurf = c0 + fcondtop = c0 + fcondbot = c0 + fbot = c0 + Tbot = c0 + Tsnice = c0 + congel = c0 + frazil = c0 + snoice = c0 + meltt = c0 + melts = c0 + meltb = c0 + meltl = c0 + dsnow = c0 + daidtt = c0 + dvidtt = c0 + dvsdtt = c0 + dagedtt = c0 + mlt_onset = c0 + frz_onset = c0 + frazil_diag= c0 + fresh_ai = c0 + fsalt_ai = c0 + fhocn_ai = c0 + fswthru_ai = c0 + fresh_da = c0 + fsalt_da = c0 + uatmT = c0 + vatmT = c0 + wlat = c0 + fsw = c0 + coszen = c0 + rdg_conv = c0 + rdg_shear = c0 + rsiden = c0 + dardg1ndt = c0 + dardg2ndt = c0 + dvirdgndt = c0 + aparticn = c0 + krdgn = c0 + ardgn = c0 + vrdgn = c0 + araftn = c0 + vraftn = c0 + aredistn = c0 + vredistn = c0 + fsurfn_f = c0 + fcondtopn_f= c0 + fsensn_f = c0 + flatn_f = c0 + evapn_f = c0 + dflatndTsfc_f = c0 + dfsurfndTsfc_f= c0 + meltsn = c0 + melttn = c0 + meltbn = c0 + congeln = c0 + snoicen = c0 + keffn_top = c0 + fsurfn = c0 + fcondtopn = c0 + fcondbotn = c0 + fsensn = c0 + flatn = c0 + albcnt = c0 + snwcnt = c0 + salinz = c0 + Tmltz = c0 + + if (grid_ice == "CD" .or. grid_ice == "C") then allocate( & - taubxN (nx_block,ny_block,max_blocks), & ! seabed stress (x) at N points (N/m^2) - taubyN (nx_block,ny_block,max_blocks), & ! seabed stress (y) at N points (N/m^2) - strairxN (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction at N points - strairyN (nx_block,ny_block,max_blocks), & ! stress on ice by air, y-direction at N points - strocnxN (nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction at N points - strocnyN (nx_block,ny_block,max_blocks), & ! ice-ocean stress, y-direction at N points - strtltxN (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, x-direction at N points - strtltyN (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, y-direction at N points - strintxN (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, x at N points (N/m^2) - strintyN (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, y at N points (N/m^2) - fmN (nx_block,ny_block,max_blocks), & ! Coriolis param. * mass in N-cell (kg/s) - TbN (nx_block,ny_block,max_blocks), & ! factor for seabed stress (landfast ice) - taubxE (nx_block,ny_block,max_blocks), & ! seabed stress (x) at E points (N/m^2) - taubyE (nx_block,ny_block,max_blocks), & ! seabed stress (y) at E points (N/m^2) - strairxE (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction at E points - strairyE (nx_block,ny_block,max_blocks), & ! stress on ice by air, y-direction at E points - strocnxE (nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction at E points - strocnyE (nx_block,ny_block,max_blocks), & ! ice-ocean stress, y-direction at E points - strtltxE (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, x-direction at E points - strtltyE (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, y-direction at E points - strintxE (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, x at E points (N/m^2) - strintyE (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, y at E points (N/m^2) - fmE (nx_block,ny_block,max_blocks), & ! Coriolis param. * mass in E-cell (kg/s) - TbE (nx_block,ny_block,max_blocks), & ! factor for seabed stress (landfast ice) - stresspT (nx_block,ny_block,max_blocks), & ! sigma11+sigma22 - stressmT (nx_block,ny_block,max_blocks), & ! sigma11-sigma22 - stress12T (nx_block,ny_block,max_blocks), & ! sigma12 - stresspU (nx_block,ny_block,max_blocks), & ! sigma11+sigma22 - stressmU (nx_block,ny_block,max_blocks), & ! sigma11-sigma22 - stress12U (nx_block,ny_block,max_blocks), & ! sigma12 - stat=ierr) - if (ierr/=0) call abort_ice('(alloc_flux): Out of memory (C or CD grid)') + taubxN (nx_block,ny_block,max_blocks), & ! seabed stress (x) at N points (N/m^2) + taubyN (nx_block,ny_block,max_blocks), & ! seabed stress (y) at N points (N/m^2) + strairxN (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction at N points + strairyN (nx_block,ny_block,max_blocks), & ! stress on ice by air, y-direction at N points + strocnxN (nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction at N points + strocnyN (nx_block,ny_block,max_blocks), & ! ice-ocean stress, y-direction at N points + strtltxN (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, x-direction at N points + strtltyN (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, y-direction at N points + strintxN (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, x at N points (N/m^2) + strintyN (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, y at N points (N/m^2) + fmN (nx_block,ny_block,max_blocks), & ! Coriolis param. * mass in N-cell (kg/s) + TbN (nx_block,ny_block,max_blocks), & ! factor for seabed stress (landfast ice) + taubxE (nx_block,ny_block,max_blocks), & ! seabed stress (x) at E points (N/m^2) + taubyE (nx_block,ny_block,max_blocks), & ! seabed stress (y) at E points (N/m^2) + strairxE (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction at E points + strairyE (nx_block,ny_block,max_blocks), & ! stress on ice by air, y-direction at E points + strocnxE (nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction at E points + strocnyE (nx_block,ny_block,max_blocks), & ! ice-ocean stress, y-direction at E points + strtltxE (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, x-direction at E points + strtltyE (nx_block,ny_block,max_blocks), & ! stress due to sea surface slope, y-direction at E points + strintxE (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, x at E points (N/m^2) + strintyE (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, y at E points (N/m^2) + fmE (nx_block,ny_block,max_blocks), & ! Coriolis param. * mass in E-cell (kg/s) + TbE (nx_block,ny_block,max_blocks), & ! factor for seabed stress (landfast ice) + stresspT (nx_block,ny_block,max_blocks), & ! sigma11+sigma22 + stressmT (nx_block,ny_block,max_blocks), & ! sigma11-sigma22 + stress12T (nx_block,ny_block,max_blocks), & ! sigma12 + stresspU (nx_block,ny_block,max_blocks), & ! sigma11+sigma22 + stressmU (nx_block,ny_block,max_blocks), & ! sigma11-sigma22 + stress12U (nx_block,ny_block,max_blocks), & ! sigma12 + stat=ierr) + if (ierr/=0) call abort_ice('(alloc_flux): Out of memory (C or CD grid)') + + taubxN = c0 + taubyN = c0 + strairxN = c0 + strairyN = c0 + strocnxN = c0 + strocnyN = c0 + strtltxN = c0 + strtltyN = c0 + strintxN = c0 + strintyN = c0 + fmN = c0 + TbN = c0 + taubxE = c0 + taubyE = c0 + strairxE = c0 + strairyE = c0 + strocnxE = c0 + strocnyE = c0 + strtltxE = c0 + strtltyE = c0 + strintxE = c0 + strintyE = c0 + fmE = c0 + TbE = c0 + stresspT = c0 + stressmT = c0 + stress12T = c0 + stresspU = c0 + stressmU = c0 + stress12U = c0 + endif ! Pond diagnostics allocate( & @@ -677,6 +907,19 @@ subroutine alloc_flux stat=ierr) if (ierr/=0) call abort_ice('(alloc_flux): Out of memory (ponds)') + dpnd_flush = c0 + dpnd_expon = c0 + dpnd_freebd = c0 + dpnd_initial = c0 + dpnd_dlid = c0 + dpnd_melt = c0 + dpnd_ridge = c0 + dpnd_flushn = c0 + dpnd_exponn = c0 + dpnd_freebdn = c0 + dpnd_initialn= c0 + dpnd_dlidn = c0 + end subroutine alloc_flux !======================================================================= diff --git a/cicecore/cicedyn/general/ice_flux_bgc.F90 b/cicecore/cicedyn/general/ice_flux_bgc.F90 index 9c07971ff..7aaaf2baa 100644 --- a/cicecore/cicedyn/general/ice_flux_bgc.F90 +++ b/cicecore/cicedyn/general/ice_flux_bgc.F90 @@ -7,6 +7,7 @@ module ice_flux_bgc use ice_kinds_mod + use ice_constants, only: c0 use ice_blocks, only: nx_block, ny_block use ice_domain_size, only: max_blocks, ncat use ice_fileunits, only: nu_diag @@ -161,6 +162,48 @@ subroutine alloc_flux_bgc stat=ierr) if (ierr/=0) call abort_ice('(alloc_flux_bgc): Out of memory') + nit = c0 + amm = c0 + sil = c0 + dmsp = c0 + dms = c0 + hum = c0 + fnit = c0 + famm = c0 + fsil = c0 + fdmsp = c0 + fdms = c0 + fhum = c0 + fdust = c0 + hin_old = c0 + dsnown = c0 + HDO_ocn = c0 + H2_16O_ocn = c0 + H2_18O_ocn = c0 + Qa_iso = c0 + Qref_iso = c0 + fiso_atm = c0 + fiso_evap = c0 + fiso_ocn = c0 + faero_atm = c0 + faero_ocn = c0 + zaeros = c0 + flux_bio_atm= c0 + flux_bio = c0 + flux_bio_ai = c0 + algalN = c0 + falgalN = c0 + doc = c0 + fdoc = c0 + don = c0 + fdon = c0 + dic = c0 + fdic = c0 + fed = c0 + fep = c0 + ffed = c0 + ffep = c0 + end subroutine alloc_flux_bgc !======================================================================= diff --git a/cicecore/cicedyn/general/ice_forcing.F90 b/cicecore/cicedyn/general/ice_forcing.F90 index e0f1b736a..d165b612a 100755 --- a/cicecore/cicedyn/general/ice_forcing.F90 +++ b/cicecore/cicedyn/general/ice_forcing.F90 @@ -209,37 +209,59 @@ subroutine alloc_forcing if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start' allocate ( & - cldf(nx_block,ny_block, max_blocks), & ! cloud fraction - fsw_data(nx_block,ny_block,2,max_blocks), & ! field values at 2 temporal data points - cldf_data(nx_block,ny_block,2,max_blocks), & - fsnow_data(nx_block,ny_block,2,max_blocks), & - Tair_data(nx_block,ny_block,2,max_blocks), & - uatm_data(nx_block,ny_block,2,max_blocks), & - vatm_data(nx_block,ny_block,2,max_blocks), & - wind_data(nx_block,ny_block,2,max_blocks), & - strax_data(nx_block,ny_block,2,max_blocks), & - stray_data(nx_block,ny_block,2,max_blocks), & - Qa_data(nx_block,ny_block,2,max_blocks), & - rhoa_data(nx_block,ny_block,2,max_blocks), & - flw_data(nx_block,ny_block,2,max_blocks), & - sst_data(nx_block,ny_block,2,max_blocks), & - sss_data(nx_block,ny_block,2,max_blocks), & - uocn_data(nx_block,ny_block,2,max_blocks), & - vocn_data(nx_block,ny_block,2,max_blocks), & - sublim_data(nx_block,ny_block,2,max_blocks), & - frain_data(nx_block,ny_block,2,max_blocks), & - topmelt_data(nx_block,ny_block,2,max_blocks,ncat), & - botmelt_data(nx_block,ny_block,2,max_blocks,ncat), & - ocn_frc_m(nx_block,ny_block, max_blocks,nfld,12), & ! ocn data for 12 months - topmelt_file(ncat), & - botmelt_file(ncat), & - wave_spectrum_data(nx_block,ny_block,nfreq,2,max_blocks), & + cldf (nx_block,ny_block, max_blocks), & ! cloud fraction + fsw_data (nx_block,ny_block,2,max_blocks), & ! field values at 2 temporal data points + cldf_data (nx_block,ny_block,2,max_blocks), & + fsnow_data (nx_block,ny_block,2,max_blocks), & + Tair_data (nx_block,ny_block,2,max_blocks), & + uatm_data (nx_block,ny_block,2,max_blocks), & + vatm_data (nx_block,ny_block,2,max_blocks), & + wind_data (nx_block,ny_block,2,max_blocks), & + strax_data (nx_block,ny_block,2,max_blocks), & + stray_data (nx_block,ny_block,2,max_blocks), & + Qa_data (nx_block,ny_block,2,max_blocks), & + rhoa_data (nx_block,ny_block,2,max_blocks), & + flw_data (nx_block,ny_block,2,max_blocks), & + sst_data (nx_block,ny_block,2,max_blocks), & + sss_data (nx_block,ny_block,2,max_blocks), & + uocn_data (nx_block,ny_block,2,max_blocks), & + vocn_data (nx_block,ny_block,2,max_blocks), & + sublim_data (nx_block,ny_block,2,max_blocks), & + frain_data (nx_block,ny_block,2,max_blocks), & + topmelt_data(nx_block,ny_block,2,max_blocks,ncat), & + botmelt_data(nx_block,ny_block,2,max_blocks,ncat), & + ocn_frc_m (nx_block,ny_block, max_blocks,nfld,12), & ! ocn data for 12 months + topmelt_file(ncat), & + botmelt_file(ncat), & + wave_spectrum_data(nx_block,ny_block,nfreq,2,max_blocks), & stat=ierr) if (ierr/=0) call abort_ice('(alloc_forcing): Out of Memory') -! initialize this, not set in box2001 (and some other forcings?) - - cldf = c0 + cldf = c0 + fsw_data = c0 + cldf_data = c0 + fsnow_data = c0 + Tair_data = c0 + uatm_data = c0 + vatm_data = c0 + wind_data = c0 + strax_data = c0 + stray_data = c0 + Qa_data = c0 + rhoa_data = c0 + flw_data = c0 + sst_data = c0 + sss_data = c0 + uocn_data = c0 + vocn_data = c0 + sublim_data = c0 + frain_data = c0 + topmelt_data = c0 + botmelt_data = c0 + ocn_frc_m = c0 + topmelt_file = '' + botmelt_file = '' + wave_spectrum_data = c0 end subroutine alloc_forcing @@ -711,13 +733,13 @@ subroutine get_forcing_atmo call ice_timer_start(timer_bound) call ice_HaloUpdate (swvdr, halo_info, & - field_loc_center, field_type_scalar) + field_loc_center, field_type_scalar, fillvalue=c0) call ice_HaloUpdate (swvdf, halo_info, & - field_loc_center, field_type_scalar) + field_loc_center, field_type_scalar, fillvalue=c0) call ice_HaloUpdate (swidr, halo_info, & - field_loc_center, field_type_scalar) + field_loc_center, field_type_scalar, fillvalue=c0) call ice_HaloUpdate (swidf, halo_info, & - field_loc_center, field_type_scalar) + field_loc_center, field_type_scalar, fillvalue=c0) call ice_timer_stop(timer_bound) call ice_timer_stop(timer_forcing) diff --git a/cicecore/cicedyn/general/ice_forcing_bgc.F90 b/cicecore/cicedyn/general/ice_forcing_bgc.F90 index 69c3ea311..d12df6417 100644 --- a/cicecore/cicedyn/general/ice_forcing_bgc.F90 +++ b/cicecore/cicedyn/general/ice_forcing_bgc.F90 @@ -65,6 +65,11 @@ subroutine alloc_forcing_bgc stat=ierr) if (ierr/=0) call abort_ice('(alloc_forcing_bgc): Out of memory') + nitdat = c0 + sildat = c0 + nit_data= c0 + sil_data= c0 + end subroutine alloc_forcing_bgc !======================================================================= diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index 92580a512..f0a6d48f1 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/cicecore/cicedyn/general/ice_init.F90 @@ -3068,14 +3068,14 @@ subroutine init_state ! Halo update on North, East faces call ice_HaloUpdate(uvelN, halo_info, & - field_loc_Nface, field_type_scalar) + field_loc_Nface, field_type_scalar, fillvalue=c0) call ice_HaloUpdate(vvelN, halo_info, & - field_loc_Nface, field_type_scalar) + field_loc_Nface, field_type_scalar, fillvalue=c0) call ice_HaloUpdate(uvelE, halo_info, & - field_loc_Eface, field_type_scalar) + field_loc_Eface, field_type_scalar, fillvalue=c0) call ice_HaloUpdate(vvelE, halo_info, & - field_loc_Eface, field_type_scalar) + field_loc_Eface, field_type_scalar, fillvalue=c0) endif diff --git a/cicecore/cicedyn/general/ice_state.F90 b/cicecore/cicedyn/general/ice_state.F90 index 82b03f2cb..92066d038 100644 --- a/cicecore/cicedyn/general/ice_state.F90 +++ b/cicecore/cicedyn/general/ice_state.F90 @@ -172,6 +172,32 @@ subroutine alloc_state stat=ierr) if (ierr/=0) call abort_ice('(alloc_state): Out of memory1') + aice = c0 + aiU = c0 + vice = c0 + vsno = c0 + aice0 = c0 + uvel = c0 + vvel = c0 + uvelE = c0 + vvelE = c0 + uvelN = c0 + vvelN = c0 + divu = c0 + shear = c0 + vort = c0 + strength = c0 + aice_init = c0 + aicen = c0 + vicen = c0 + vsnon = c0 + aicen_init = c0 + vicen_init = c0 + vsnon_init = c0 + Tsfcn_init = c0 + trcr = c0 + trcrn = c0 + allocate ( & trcr_depend(ntrcr) , & ! n_trcr_strata(ntrcr) , & ! number of underlying tracer layers @@ -184,12 +210,6 @@ subroutine alloc_state n_trcr_strata = 0 nt_strata = 0 trcr_base = c0 - aicen = c0 - aicen_init = c0 - vicen = c0 - vicen_init = c0 - vsnon = c0 - vsnon_init = c0 end subroutine alloc_state diff --git a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 index 11cd0d2e1..5a690d490 100644 --- a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 @@ -64,7 +64,7 @@ module ice_boundary use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use ice_blocks, only: nx_block, ny_block, nghost, & - nblocks_tot, ice_blocksNorth, & + nblocks_tot, nblocks_x, nblocks_y, ice_blocksNorth, & ice_blocksSouth, ice_blocksEast, ice_blocksWest, & ice_blocksEast2, ice_blocksWest2, & ice_blocksNorthEast, ice_blocksNorthWest, & @@ -106,6 +106,10 @@ module ice_boundary sendAddr, &! src addresses for each sent message recvAddr ! dst addresses for each recvd message + character (char_len) :: & + nsBoundaryType, &! type of boundary to use in logical ns dir + ewBoundaryType ! type of boundary to use in logical ew dir + end type public :: ice_HaloCreate, & @@ -252,6 +256,8 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & if (my_task >= numProcs) return halo%communicator = communicator + halo%ewBoundaryType = ewBoundaryType + halo%nsBoundaryType = nsBoundaryType blockSizeX = nx_block - 2*nghost blockSizeY = ny_block - 2*nghost @@ -1239,6 +1245,7 @@ subroutine ice_HaloUpdate2DR8(array, halo, & integer (int_kind) :: & i,j,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) iSrc,jSrc, &! source addresses for message @@ -1261,6 +1268,8 @@ subroutine ice_HaloUpdate2DR8(array, halo, & x1,x2,xavg ! scalars for enforcing symmetry at U pts logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter, &! fill outer boundary ns ltripoleOnly ! local tripoleOnly value integer (int_kind) :: len ! length of messages @@ -1290,8 +1299,20 @@ subroutine ice_HaloUpdate2DR8(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_dbl_kind endif @@ -1367,29 +1388,77 @@ subroutine ice_HaloUpdate2DR8(array, halo, & !----------------------------------------------------------------------- ! -! while messages are being communicated, fill out halo region +! While messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- - if (ltripoleOnly) then - ! skip fill, not needed since tripole seam always exists if running - ! on tripole grid and set tripoleOnly flag - else + if (.not. ltripoleOnly) then + ! tripoleOnly skip fill, do not overwrite any values in interior as they may + ! already be set and filling tripole is not necessary + + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,iblk) = fill - array(1:nx_block, jhi+j,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,iblk) = fill - array(ihi+i, 1:ny_block,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, iblk) = fill + enddo + enddo + endif + enddo ! iblk endif !----------------------------------------------------------------------- @@ -1683,6 +1752,7 @@ subroutine ice_HaloUpdate2DR4(array, halo, & integer (int_kind) :: & i,j,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) iSrc,jSrc, &! source addresses for message @@ -1704,6 +1774,10 @@ subroutine ice_HaloUpdate2DR4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + integer (int_kind) :: len ! length of messages character(len=*), parameter :: subname = '(ice_HaloUpdate2DR4)' @@ -1731,8 +1805,20 @@ subroutine ice_HaloUpdate2DR4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_real_kind endif @@ -1804,23 +1890,71 @@ subroutine ice_HaloUpdate2DR4(array, halo, & ! ! while messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,iblk) = fill - array(1:nx_block, jhi+j,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,iblk) = fill - array(ihi+i, 1:ny_block,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -2098,6 +2232,7 @@ subroutine ice_HaloUpdate2DI4(array, halo, & integer (int_kind) :: & i,j,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) iSrc,jSrc, &! source addresses for message @@ -2119,6 +2254,10 @@ subroutine ice_HaloUpdate2DI4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + integer (int_kind) :: len ! length of messages character(len=*), parameter :: subname = '(ice_HaloUpdate2DI4)' @@ -2146,8 +2285,20 @@ subroutine ice_HaloUpdate2DI4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0_int_kind endif @@ -2219,23 +2370,71 @@ subroutine ice_HaloUpdate2DI4(array, halo, & ! ! while messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,iblk) = fill - array(1:nx_block, jhi+j,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,iblk) = fill - array(ihi+i, 1:ny_block,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -2593,6 +2792,7 @@ subroutine ice_HaloUpdate3DR8(array, halo, & integer (int_kind) :: & i,j,k,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) nz, &! size of array in 3rd dimension @@ -2615,6 +2815,10 @@ subroutine ice_HaloUpdate3DR8(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + real (dbl_kind), dimension(:,:), allocatable :: & bufSend, bufRecv ! 3d send,recv buffers @@ -2648,8 +2852,20 @@ subroutine ice_HaloUpdate3DR8(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_dbl_kind endif @@ -2742,23 +2958,71 @@ subroutine ice_HaloUpdate3DR8(array, halo, & ! ! while messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,iblk) = fill - array(1:nx_block, jhi+j,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,iblk) = fill - array(ihi+i, 1:ny_block,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -3067,6 +3331,7 @@ subroutine ice_HaloUpdate3DR4(array, halo, & integer (int_kind) :: & i,j,k,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) nz, &! size of array in 3rd dimension @@ -3089,6 +3354,10 @@ subroutine ice_HaloUpdate3DR4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + real (real_kind), dimension(:,:), allocatable :: & bufSend, bufRecv ! 3d send,recv buffers @@ -3122,8 +3391,20 @@ subroutine ice_HaloUpdate3DR4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_real_kind endif @@ -3216,23 +3497,71 @@ subroutine ice_HaloUpdate3DR4(array, halo, & ! ! while messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,iblk) = fill - array(1:nx_block, jhi+j,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,iblk) = fill - array(ihi+i, 1:ny_block,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -3541,6 +3870,7 @@ subroutine ice_HaloUpdate3DI4(array, halo, & integer (int_kind) :: & i,j,k,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) nz, &! size of array in 3rd dimension @@ -3563,6 +3893,10 @@ subroutine ice_HaloUpdate3DI4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + integer (int_kind), dimension(:,:), allocatable :: & bufSend, bufRecv ! 3d send,recv buffers @@ -3596,8 +3930,20 @@ subroutine ice_HaloUpdate3DI4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0_int_kind endif @@ -3690,23 +4036,71 @@ subroutine ice_HaloUpdate3DI4(array, halo, & ! ! while messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,iblk) = fill - array(1:nx_block, jhi+j,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,iblk) = fill - array(ihi+i, 1:ny_block,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -4015,6 +4409,7 @@ subroutine ice_HaloUpdate4DR8(array, halo, & integer (int_kind) :: & i,j,k,l,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) nz, nt, &! size of array in 3rd,4th dimensions @@ -4037,6 +4432,10 @@ subroutine ice_HaloUpdate4DR8(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + real (dbl_kind), dimension(:,:), allocatable :: & bufSend, bufRecv ! 4d send,recv buffers @@ -4070,8 +4469,20 @@ subroutine ice_HaloUpdate4DR8(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_dbl_kind endif @@ -4168,23 +4579,71 @@ subroutine ice_HaloUpdate4DR8(array, halo, & ! ! while messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,:,iblk) = fill - array(1:nx_block, jhi+j,:,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,:,iblk) = fill - array(ihi+i, 1:ny_block,:,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -4513,6 +4972,7 @@ subroutine ice_HaloUpdate4DR4(array, halo, & integer (int_kind) :: & i,j,k,l,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) nz, nt, &! size of array in 3rd,4th dimensions @@ -4535,6 +4995,10 @@ subroutine ice_HaloUpdate4DR4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + real (real_kind), dimension(:,:), allocatable :: & bufSend, bufRecv ! 4d send,recv buffers @@ -4568,8 +5032,20 @@ subroutine ice_HaloUpdate4DR4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_real_kind endif @@ -4666,23 +5142,71 @@ subroutine ice_HaloUpdate4DR4(array, halo, & ! ! while messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,:,iblk) = fill - array(1:nx_block, jhi+j,:,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,:,iblk) = fill - array(ihi+i, 1:ny_block,:,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -5011,6 +5535,7 @@ subroutine ice_HaloUpdate4DI4(array, halo, & integer (int_kind) :: & i,j,k,l,n,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices ierr, &! error or status flag for MPI,alloc nxGlobal, &! global domain size in x (tripole) nz, nt, &! size of array in 3rd,4th dimensions @@ -5033,6 +5558,10 @@ subroutine ice_HaloUpdate4DI4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + integer (int_kind), dimension(:,:), allocatable :: & bufSend, bufRecv ! 4d send,recv buffers @@ -5066,8 +5595,20 @@ subroutine ice_HaloUpdate4DI4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0_int_kind endif @@ -5164,23 +5705,71 @@ subroutine ice_HaloUpdate4DI4(array, halo, & ! ! while messages are being communicated, fill out halo region ! needed for masked halos to ensure halo values are filled for -! halo grid cells that are not updated +! halo grid cells that are not updated except in cases where +! you don't want to overwrite those halos ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,:,iblk) = fill - array(1:nx_block, jhi+j,:,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,:,iblk) = fill - array(ihi+i, 1:ny_block,:,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -7027,15 +7616,15 @@ end subroutine ice_HaloMsgCreate subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) -! This subroutine extrapolates ARRAY values into the first row or column -! of ghost cells, and is intended for grid variables whose ghost cells +! This subroutine extrapolates ARRAY values into the ghost cells, +! and is intended for grid variables whose ghost cells ! would otherwise be set using the default boundary conditions (Dirichlet ! or Neumann). -! Note: This routine will need to be modified for nghost > 1. -! We assume padding occurs only on east and north edges. ! ! This is the specific interface for double precision arrays ! corresponding to the generic interface ice_HaloExtrapolate +! +! T.Craig, Oct 2025 - extend to nghost > 1 use ice_blocks, only: block, nblocks_x, nblocks_y, get_block use ice_constants, only: c2 @@ -7058,8 +7647,9 @@ subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) !----------------------------------------------------------------------- integer (int_kind) :: & - i,j,iblk, &! dummy loop indices - numBlocks, &! number of local blocks + i,j,n,iblk,ii,jj, &! dummy loop indices + ilo,ihi,jlo,jhi, &! active block indices + numBlocks, &! number of local blocks blockID, &! block location ibc ! ghost cell column or row @@ -7067,6 +7657,7 @@ subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) this_block ! block info for current block character(len=*), parameter :: subname = '(ice_HaloExtrapolate2DR8)' + !----------------------------------------------------------------------- ! ! Linear extrapolation @@ -7079,32 +7670,40 @@ subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) do iblk = 1, numBlocks call ice_distributionGetBlockID(dist, iblk, blockID) this_block = get_block(blockID, blockID) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi if (this_block%iblock == 1) then ! west edge if (trim(ew_bndy_type) /= 'cyclic') then + do n = 1, nghost + ii = ilo - n ! gridcell to extrapolate to do j = 1, ny_block - ARRAY(1,j,iblk) = c2*ARRAY(2,j,iblk) - ARRAY(3,j,iblk) + ARRAY(ii,j,iblk) = c2*ARRAY(ii+1,j,iblk) - ARRAY(ii+2,j,iblk) + enddo enddo endif endif if (this_block%iblock == nblocks_x) then ! east edge if (trim(ew_bndy_type) /= 'cyclic') then - ! locate ghost cell column (avoid padding) - ibc = nx_block - do i = nx_block, nghost + 1, -1 - if (this_block%i_glob(i) == 0) ibc = ibc - 1 - enddo + do n = 1, nghost + ii = ihi + n ! gridcell to extrapolate to do j = 1, ny_block - ARRAY(ibc,j,iblk) = c2*ARRAY(ibc-1,j,iblk) - ARRAY(ibc-2,j,iblk) + ARRAY(ii,j,iblk) = c2*ARRAY(ii-1,j,iblk) - ARRAY(ii-2,j,iblk) + enddo enddo endif endif if (this_block%jblock == 1) then ! south edge if (trim(ns_bndy_type) /= 'cyclic') then + do n = 1, nghost + jj = jlo - n ! gridcell to extrapolate to do i = 1, nx_block - ARRAY(i,1,iblk) = c2*ARRAY(i,2,iblk) - ARRAY(i,3,iblk) + ARRAY(i,jj,iblk) = c2*ARRAY(i,jj+1,iblk) - ARRAY(i,jj+2,iblk) + enddo enddo endif endif @@ -7113,13 +7712,11 @@ subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) if (trim(ns_bndy_type) /= 'cyclic' .and. & trim(ns_bndy_type) /= 'tripole' .and. & trim(ns_bndy_type) /= 'tripoleT' ) then - ! locate ghost cell column (avoid padding) - ibc = ny_block - do j = ny_block, nghost + 1, -1 - if (this_block%j_glob(j) == 0) ibc = ibc - 1 - enddo + do n = 1, nghost + jj = jhi + n ! gridcell to extrapolate to do i = 1, nx_block - ARRAY(i,ibc,iblk) = c2*ARRAY(i,ibc-1,iblk) - ARRAY(i,ibc-2,iblk) + ARRAY(i,jj,iblk) = c2*ARRAY(i,jj-1,iblk) - ARRAY(i,jj-2,iblk) + enddo enddo endif endif diff --git a/cicecore/cicedyn/infrastructure/comm/mpi/ice_gather_scatter.F90 b/cicecore/cicedyn/infrastructure/comm/mpi/ice_gather_scatter.F90 index cfb98befe..cd3b4afba 100644 --- a/cicecore/cicedyn/infrastructure/comm/mpi/ice_gather_scatter.F90 +++ b/cicecore/cicedyn/infrastructure/comm/mpi/ice_gather_scatter.F90 @@ -1643,74 +1643,35 @@ subroutine scatter_global_dbl(ARRAY, ARRAY_G, src_task, dst_dist, & msg_buffer = c0 this_block = get_block(n,n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - msg_buffer(i,j-yoffset2) = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + msg_buffer(i,j-yoffset2) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + msg_buffer(i,j) = ARRAY_G(isrc,jsrc) + end do + endif + end do call MPI_SEND(msg_buffer, nx_block*ny_block, & mpiR8, dst_dist%blockLocation(n)-1, 3*mpitag_gs+n, & @@ -1728,75 +1689,35 @@ subroutine scatter_global_dbl(ARRAY, ARRAY_G, src_task, dst_dist, & dst_block = dst_dist%blockLocalID(n) this_block = get_block(n,n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - ARRAY(i,j-yoffset2,dst_block) & - = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + ARRAY(i,j-yoffset2,dst_block) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + ARRAY(i,j,dst_block) = ARRAY_G(isrc,jsrc) + end do + endif + end do endif end do @@ -1832,7 +1753,7 @@ subroutine scatter_global_dbl(ARRAY, ARRAY_G, src_task, dst_dist, & endif !----------------------------------------------------------------- - ! Ensure unused ghost cell values are 0 + ! Ensure/reset ghost cell values are 0 for noupdate !----------------------------------------------------------------- if (field_loc == field_loc_noupdate) then @@ -2029,74 +1950,35 @@ subroutine scatter_global_real(ARRAY, ARRAY_G, src_task, dst_dist, & msg_buffer = 0._real_kind this_block = get_block(n,n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - msg_buffer(i,j-yoffset2) = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + msg_buffer(i,j-yoffset2) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + msg_buffer(i,j) = ARRAY_G(isrc,jsrc) + end do + endif + end do call MPI_SEND(msg_buffer, nx_block*ny_block, & mpiR4, dst_dist%blockLocation(n)-1, 3*mpitag_gs+n, & @@ -2114,75 +1996,35 @@ subroutine scatter_global_real(ARRAY, ARRAY_G, src_task, dst_dist, & dst_block = dst_dist%blockLocalID(n) this_block = get_block(n,n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - ARRAY(i,j-yoffset2,dst_block) & - = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + ARRAY(i,j-yoffset2,dst_block) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + ARRAY(i,j,dst_block) = ARRAY_G(isrc,jsrc) + end do + endif + end do endif end do @@ -2218,7 +2060,7 @@ subroutine scatter_global_real(ARRAY, ARRAY_G, src_task, dst_dist, & endif !----------------------------------------------------------------- - ! Ensure unused ghost cell values are 0 + ! Ensure/reset ghost cell values are 0 for noupdate !----------------------------------------------------------------- if (field_loc == field_loc_noupdate) then @@ -2415,74 +2257,35 @@ subroutine scatter_global_int(ARRAY, ARRAY_G, src_task, dst_dist, & msg_buffer = 0 this_block = get_block(n,n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - msg_buffer(i,j) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - msg_buffer(i,j-yoffset2) = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + msg_buffer(i,j-yoffset2) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + msg_buffer(i,j) = ARRAY_G(isrc,jsrc) + end do + endif + end do call MPI_SEND(msg_buffer, nx_block*ny_block, & mpi_integer, dst_dist%blockLocation(n)-1, 3*mpitag_gs+n, & @@ -2500,75 +2303,35 @@ subroutine scatter_global_int(ARRAY, ARRAY_G, src_task, dst_dist, & dst_block = dst_dist%blockLocalID(n) this_block = get_block(n,n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - ARRAY(i,j-yoffset2,dst_block) & - = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + ARRAY(i,j-yoffset2,dst_block) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + ARRAY(i,j,dst_block) = ARRAY_G(isrc,jsrc) + end do + endif + end do endif end do @@ -2604,7 +2367,7 @@ subroutine scatter_global_int(ARRAY, ARRAY_G, src_task, dst_dist, & endif !----------------------------------------------------------------- - ! Ensure unused ghost cell values are 0 + ! Ensure/reset ghost cell values are 0 for noupdate !----------------------------------------------------------------- if (field_loc == field_loc_noupdate) then @@ -2681,7 +2444,7 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) integer (int_kind) :: & i,j,n, &! dummy loop indices iblk, jblk, &! block indices - iglb, jglb, &! global indices + isrc, jsrc, &! global indices nrecvs, &! actual number of messages received dst_block, &! location of block in dst array ierr ! MPI error flag @@ -2748,13 +2511,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! southwest corner iblk = i jblk = j - iglb = this_block%i_glob(this_block%ilo)+i-1 - jglb = j - msg_buffer(iblk,jblk) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ilo)+i-1 + jsrc = j + msg_buffer(iblk,jblk) = ARRAY_G(isrc,jsrc) ! southeast corner iblk = this_block%ihi+i - iglb = this_block%i_glob(this_block%ihi)+nghost+i - msg_buffer(iblk,jblk) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ihi)+nghost+i + msg_buffer(iblk,jblk) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -2769,13 +2532,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northwest corner iblk = i jblk = this_block%jhi+j - iglb = this_block%i_glob(this_block%ilo)+i-1 - jglb = ny_global+nghost+j - msg_buffer(iblk,jblk) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ilo)+i-1 + jsrc = ny_global+nghost+j + msg_buffer(iblk,jblk) = ARRAY_G(isrc,jsrc) ! northeast corner iblk = this_block%ihi+i - iglb = this_block%i_glob(this_block%ihi)+nghost+i - msg_buffer(iblk,jblk) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ihi)+nghost+i + msg_buffer(iblk,jblk) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -2791,13 +2554,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northwest corner iblk = i jblk = this_block%jhi+j - iglb = i - jglb = this_block%j_glob(this_block%jhi)+nghost+j - msg_buffer(iblk,jblk) = ARRAY_G(iglb,jglb) + isrc = i + jsrc = this_block%j_glob(this_block%jhi)+nghost+j + msg_buffer(iblk,jblk) = ARRAY_G(isrc,jsrc) ! southwest corner jblk = j - jglb = this_block%j_glob(this_block%jlo)+j-1 - msg_buffer(iblk,jblk) = ARRAY_G(iglb,jglb) + jsrc = this_block%j_glob(this_block%jlo)+j-1 + msg_buffer(iblk,jblk) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -2814,13 +2577,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northeast corner iblk = this_block%ihi+i jblk = this_block%jhi+j - iglb = nx_global+nghost+i - jglb = this_block%j_glob(this_block%jhi)+nghost+j - msg_buffer(iblk,jblk) = ARRAY_G(iglb,jglb) + isrc = nx_global+nghost+i + jsrc = this_block%j_glob(this_block%jhi)+nghost+j + msg_buffer(iblk,jblk) = ARRAY_G(isrc,jsrc) ! southeast corner jblk = j - jglb = this_block%j_glob(this_block%jlo)+j-1 - msg_buffer(iblk,jblk) = ARRAY_G(iglb,jglb) + jsrc = this_block%j_glob(this_block%jlo)+j-1 + msg_buffer(iblk,jblk) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -2861,13 +2624,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! southwest corner iblk = i jblk = j - iglb = this_block%i_glob(this_block%ilo)+i-1 - jglb = j - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ilo)+i-1 + jsrc = j + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) ! southeast corner iblk = this_block%ihi+i - iglb = this_block%i_glob(this_block%ihi)+nghost+i - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ihi)+nghost+i + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -2882,13 +2645,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northwest corner iblk = i jblk = this_block%jhi+j - iglb = this_block%i_glob(this_block%ilo)+i-1 - jglb = ny_global+nghost+j - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ilo)+i-1 + jsrc = ny_global+nghost+j + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) ! northeast corner iblk = this_block%ihi+i - iglb = this_block%i_glob(this_block%ihi)+nghost+i - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ihi)+nghost+i + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -2904,13 +2667,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northwest corner iblk = i jblk = this_block%jhi+j - iglb = i - jglb = this_block%j_glob(this_block%jhi)+nghost+j - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = i + jsrc = this_block%j_glob(this_block%jhi)+nghost+j + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) ! southwest corner jblk = j - jglb = this_block%j_glob(this_block%jlo)+j-1 - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + jsrc = this_block%j_glob(this_block%jlo)+j-1 + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -2927,17 +2690,16 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northeast corner iblk = this_block%ihi+i jblk = this_block%jhi+j - iglb = nx_global+nghost+i - jglb = this_block%j_glob(this_block%jhi)+nghost+j - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = nx_global+nghost+i + jsrc = this_block%j_glob(this_block%jhi)+nghost+j + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) ! southeast corner jblk = j - jglb = this_block%j_glob(this_block%jlo)+j-1 - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + jsrc = this_block%j_glob(this_block%jlo)+j-1 + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) enddo enddo endif - endif end do @@ -3071,70 +2833,30 @@ subroutine scatter_global_stress(ARRAY, ARRAY_G1, ARRAY_G2, & msg_buffer = c0 this_block = get_block(n,n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - msg_buffer(i,j) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - msg_buffer(i,j) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - msg_buffer(i,j) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - jsrc = ny_global + yoffset + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - msg_buffer(i,j) = isign * ARRAY_G2(isrc,jsrc) - endif - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + jsrc = ny_global + yoffset + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + msg_buffer(i,j) = isign * ARRAY_G2(isrc,jsrc) + endif + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + msg_buffer(i,j) = ARRAY_G1(isrc,jsrc) + end do + endif + end do call MPI_SEND(msg_buffer, nx_block*ny_block, & mpiR8, dst_dist%blockLocation(n)-1, 3*mpitag_gs+n, & @@ -3152,75 +2874,35 @@ subroutine scatter_global_stress(ARRAY, ARRAY_G1, ARRAY_G2, & dst_block = dst_dist%blockLocalID(n) this_block = get_block(n,n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - ARRAY(i,j,dst_block) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - ARRAY(i,j-yoffset2,dst_block) & - = isign * ARRAY_G2(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + ARRAY(i,j-yoffset2,dst_block) & + = isign * ARRAY_G2(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + ARRAY(i,j,dst_block) = ARRAY_G1(isrc,jsrc) + end do + endif + end do endif end do diff --git a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 index b9ac8fe33..f185da3c5 100644 --- a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 @@ -31,7 +31,7 @@ module ice_boundary use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use ice_blocks, only: nx_block, ny_block, nghost, & - nblocks_tot, ice_blocksNorth, & + nblocks_tot, ice_blocksNorth, nblocks_x, nblocks_y, & ice_blocksSouth, ice_blocksEast, ice_blocksWest, & ice_blocksEast2, ice_blocksWest2, & ice_blocksNorthEast, ice_blocksNorthWest, & @@ -61,6 +61,10 @@ module ice_boundary srcLocalAddr, &! src addresses for each local copy dstLocalAddr ! dst addresses for each local copy + character (char_len) :: & + nsBoundaryType, &! type of boundary to use in logical ns dir + ewBoundaryType ! type of boundary to use in logical ew dir + end type public :: ice_HaloCreate, & @@ -177,6 +181,8 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & if (my_task >= numProcs) return halo%communicator = communicator + halo%ewBoundaryType = ewBoundaryType + halo%nsBoundaryType = nsBoundaryType blockSizeX = nx_block - 2*nghost blockSizeY = ny_block - 2*nghost @@ -659,7 +665,7 @@ subroutine ice_HaloUpdate2DR8(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 2d horizontal arrays of double precision. type (ice_halo), intent(in) :: & @@ -690,9 +696,10 @@ subroutine ice_HaloUpdate2DR8(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) iSrc,jSrc, &! source addresses for message iDst,jDst, &! dest addresses for message @@ -706,6 +713,8 @@ subroutine ice_HaloUpdate2DR8(array, halo, & x1,x2,xavg ! scalars for enforcing symmetry at U pts logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter, &! fill outer boundary ns ltripoleOnly ! local tripoleOnly value character(len=*), parameter :: subname = '(ice_HaloUpdate2DR8)' @@ -733,8 +742,20 @@ subroutine ice_HaloUpdate2DR8(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic or tripole + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_dbl_kind endif @@ -753,29 +774,78 @@ subroutine ice_HaloUpdate2DR8(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- - if (ltripoleOnly) then - ! skip fill, not needed since tripole seam always exists if running - ! on tripole grid and set tripoleOnly flag - else + if (.not. ltripoleOnly) then + ! tripoleOnly skip fill, do not overwrite any values in interior as they may + ! already be set and filling tripole is not necessary + + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,iblk) = fill - array(1:nx_block, jhi+j,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,iblk) = fill - array(ihi+i, 1:ny_block,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, iblk) = fill + enddo + enddo + endif + enddo ! iblk endif !----------------------------------------------------------------------- @@ -994,7 +1064,7 @@ subroutine ice_HaloUpdate2DR4(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 2d horizontal arrays of single precision. type (ice_halo), intent(in) :: & @@ -1022,9 +1092,10 @@ subroutine ice_HaloUpdate2DR4(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) iSrc,jSrc, &! source addresses for message iDst,jDst, &! dest addresses for message @@ -1037,6 +1108,10 @@ subroutine ice_HaloUpdate2DR4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + character(len=*), parameter :: subname = '(ice_HaloUpdate2DR4)' !----------------------------------------------------------------------- @@ -1062,8 +1137,20 @@ subroutine ice_HaloUpdate2DR4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_real_kind endif @@ -1076,25 +1163,74 @@ subroutine ice_HaloUpdate2DR4(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,iblk) = fill - array(1:nx_block, jhi+j,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,iblk) = fill - array(ihi+i, 1:ny_block,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -1303,7 +1439,7 @@ subroutine ice_HaloUpdate2DI4(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 2d horizontal integer arrays. type (ice_halo), intent(in) :: & @@ -1331,9 +1467,10 @@ subroutine ice_HaloUpdate2DI4(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) iSrc,jSrc, &! source addresses for message iDst,jDst, &! dest addresses for message @@ -1346,6 +1483,10 @@ subroutine ice_HaloUpdate2DI4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + character(len=*), parameter :: subname = '(ice_HaloUpdate2DI4)' !----------------------------------------------------------------------- @@ -1371,8 +1512,20 @@ subroutine ice_HaloUpdate2DI4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0_int_kind endif @@ -1385,25 +1538,74 @@ subroutine ice_HaloUpdate2DI4(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,iblk) = fill - array(1:nx_block, jhi+j,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,iblk) = fill - array(ihi+i, 1:ny_block,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -1692,7 +1894,7 @@ subroutine ice_HaloUpdate3DR8(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 3d horizontal arrays of double precision. type (ice_halo), intent(in) :: & @@ -1720,9 +1922,10 @@ subroutine ice_HaloUpdate3DR8(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,k,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) nz, &! size of array in 3rd dimension iSrc,jSrc, &! source addresses for message @@ -1736,6 +1939,10 @@ subroutine ice_HaloUpdate3DR8(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + real (dbl_kind), dimension(:,:,:), allocatable :: & bufTripole ! 3d tripole buffer @@ -1764,8 +1971,20 @@ subroutine ice_HaloUpdate3DR8(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_dbl_kind endif @@ -1781,25 +2000,74 @@ subroutine ice_HaloUpdate3DR8(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,iblk) = fill - array(1:nx_block, jhi+j,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,iblk) = fill - array(ihi+i, 1:ny_block,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -2027,7 +2295,7 @@ subroutine ice_HaloUpdate3DR4(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 3d horizontal arrays of single precision. type (ice_halo), intent(in) :: & @@ -2055,9 +2323,10 @@ subroutine ice_HaloUpdate3DR4(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,k,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) nz, &! size of array in 3rd dimension iSrc,jSrc, &! source addresses for message @@ -2071,6 +2340,10 @@ subroutine ice_HaloUpdate3DR4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + real (real_kind), dimension(:,:,:), allocatable :: & bufTripole ! 3d tripole buffer @@ -2099,8 +2372,20 @@ subroutine ice_HaloUpdate3DR4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_real_kind endif @@ -2116,25 +2401,74 @@ subroutine ice_HaloUpdate3DR4(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,iblk) = fill - array(1:nx_block, jhi+j,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,iblk) = fill - array(ihi+i, 1:ny_block,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -2362,7 +2696,7 @@ subroutine ice_HaloUpdate3DI4(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 3d horizontal arrays of double precision. type (ice_halo), intent(in) :: & @@ -2390,9 +2724,10 @@ subroutine ice_HaloUpdate3DI4(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,k,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) nz, &! size of array in 3rd dimension iSrc,jSrc, &! source addresses for message @@ -2406,6 +2741,10 @@ subroutine ice_HaloUpdate3DI4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + integer (int_kind), dimension(:,:,:), allocatable :: & bufTripole ! 3d tripole buffer @@ -2434,8 +2773,20 @@ subroutine ice_HaloUpdate3DI4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0_int_kind endif @@ -2451,25 +2802,74 @@ subroutine ice_HaloUpdate3DI4(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,iblk) = fill - array(1:nx_block, jhi+j,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,iblk) = fill - array(ihi+i, 1:ny_block,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -2697,7 +3097,7 @@ subroutine ice_HaloUpdate4DR8(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 4d horizontal arrays of double precision. type (ice_halo), intent(in) :: & @@ -2725,9 +3125,10 @@ subroutine ice_HaloUpdate4DR8(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,k,l,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) nz, nt, &! size of array in 3rd,4th dimensions iSrc,jSrc, &! source addresses for message @@ -2741,6 +3142,10 @@ subroutine ice_HaloUpdate4DR8(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + real (dbl_kind), dimension(:,:,:,:), allocatable :: & bufTripole ! 4d tripole buffer @@ -2769,8 +3174,20 @@ subroutine ice_HaloUpdate4DR8(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_dbl_kind endif @@ -2787,25 +3204,74 @@ subroutine ice_HaloUpdate4DR8(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,:,iblk) = fill - array(1:nx_block, jhi+j,:,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,:,iblk) = fill - array(ihi+i, 1:ny_block,:,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -3049,7 +3515,7 @@ subroutine ice_HaloUpdate4DR4(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 4d horizontal arrays of single precision. type (ice_halo), intent(in) :: & @@ -3077,9 +3543,10 @@ subroutine ice_HaloUpdate4DR4(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,k,l,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) nz, nt, &! size of array in 3rd,4th dimensions iSrc,jSrc, &! source addresses for message @@ -3093,6 +3560,10 @@ subroutine ice_HaloUpdate4DR4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + real (real_kind), dimension(:,:,:,:), allocatable :: & bufTripole ! 4d tripole buffer @@ -3121,8 +3592,20 @@ subroutine ice_HaloUpdate4DR4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0.0_real_kind endif @@ -3139,25 +3622,74 @@ subroutine ice_HaloUpdate4DR4(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,:,iblk) = fill - array(1:nx_block, jhi+j,:,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,:,iblk) = fill - array(ihi+i, 1:ny_block,:,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -3401,7 +3933,7 @@ subroutine ice_HaloUpdate4DI4(array, halo, & ! This routine updates ghost cells for an input array and is a ! member of a group of routines under the generic interface -! POP\_HaloUpdate. This routine is the specific interface +! ice\_HaloUpdate. This routine is the specific interface ! for 4d horizontal integer arrays. type (ice_halo), intent(in) :: & @@ -3429,9 +3961,10 @@ subroutine ice_HaloUpdate4DI4(array, halo, & ! !----------------------------------------------------------------------- - integer (int_kind) :: & + integer (int_kind) :: & i,j,k,l,nmsg, &! dummy loop indices iblk,ilo,ihi,jlo,jhi, &! block sizes for fill + iblock,jblock, &! global block indices nxGlobal, &! global domain size in x (tripole) nz, nt, &! size of array in 3rd,4th dimensions iSrc,jSrc, &! source addresses for message @@ -3445,6 +3978,10 @@ subroutine ice_HaloUpdate4DI4(array, halo, & fill, &! value to use for unknown points x1,x2,xavg ! scalars for enforcing symmetry at U pts + logical (log_kind) :: & + ewfillouter, &! fill outer boundary ew + nsfillouter ! fill outer boundary ns + integer (int_kind), dimension(:,:,:,:), allocatable :: & bufTripole ! 4d tripole buffer @@ -3473,8 +4010,20 @@ subroutine ice_HaloUpdate4DI4(array, halo, & ! !----------------------------------------------------------------------- + ewfillouter = .false. + nsfillouter = .false. + + ! fill outer boundary if cyclic + if (halo%ewBoundaryType == 'cyclic') ewfillouter=.true. + if (halo%nsBoundaryType == 'tripole' .or. & + halo%nsBoundaryType == 'tripoleT' .or. & + halo%nsBoundaryType == 'cyclic') nsfillouter=.true. + if (present(fillValue)) then fill = fillValue + ! always fill outer boundary if fillValue is passed + ewfillouter = .true. + nsfillouter = .true. else fill = 0_int_kind endif @@ -3491,25 +4040,74 @@ subroutine ice_HaloUpdate4DI4(array, halo, & !----------------------------------------------------------------------- ! -! fill out halo region -! needed for masked halos to ensure halo values are filled for +! Fill out halo region +! Needed for masked halos to ensure halo values are filled for ! halo grid cells that are not updated +! In general, do NOT fill outer boundary for open boundary conditions +! because do not want to overwrite existing data ! !----------------------------------------------------------------------- + ! fill outer boundary as needed + ! only fill corners if both edges are being filled do iblk = 1, halo%numLocalBlocks - call get_block_parameter(halo%blockGlobalID(iblk), & - ilo=ilo, ihi=ihi, & - jlo=jlo, jhi=jhi) - do j = 1,nghost - array(1:nx_block, jlo-j,:,:,iblk) = fill - array(1:nx_block, jhi+j,:,:,iblk) = fill - enddo - do i = 1,nghost - array(ilo-i, 1:ny_block,:,:,iblk) = fill - array(ihi+i, 1:ny_block,:,:,iblk) = fill - enddo - enddo + call get_block_parameter(halo%blockGlobalID(iblk), & + ilo=ilo, ihi=ihi, & + jlo=jlo, jhi=jhi, & + iblock=iblock, jblock=jblock) + if (ewfillouter .or. iblock > 1) then ! west edge + do i = 1,nghost + array(ilo-i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (ewfillouter .or. iblock < nblocks_x) then ! east edge + do i = 1,nghost + array(ihi+i, jlo:jhi, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock > 1) then ! south edge + do j = 1,nghost + array(ilo:ihi, jlo-j, :, :, iblk) = fill + enddo + endif + if (nsfillouter .or. jblock < nblocks_y) then ! north edge + do j = 1,nghost + array(ilo:ihi, jhi+j, :, :, iblk) = fill + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! southwest corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock > 1) .and. & ! northwest corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ilo-i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! southeast corner + (nsfillouter .or. jblock > 1)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jlo-j, :, :, iblk) = fill + enddo + enddo + endif + if ((ewfillouter .or. iblock < nblocks_x) .and. & ! northeast corner + (nsfillouter .or. jblock < nblocks_y)) then + do j = 1,nghost + do i = 1,nghost + array(ihi+i, jhi+j, :, :, iblk) = fill + enddo + enddo + endif + enddo ! iblk !----------------------------------------------------------------------- ! @@ -4778,15 +5376,15 @@ end subroutine ice_HaloMsgCreate subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) -! This subroutine extrapolates ARRAY values into the first row or column -! of ghost cells, and is intended for grid variables whose ghost cells +! This subroutine extrapolates ARRAY values into the ghost cells, +! and is intended for grid variables whose ghost cells ! would otherwise be set using the default boundary conditions (Dirichlet ! or Neumann). -! Note: This routine will need to be modified for nghost > 1. -! We assume padding occurs only on east and north edges. ! ! This is the specific interface for double precision arrays ! corresponding to the generic interface ice_HaloExtrapolate +! +! T.Craig, Oct 2025 - extend to nghost > 1 use ice_blocks, only: block, nblocks_x, nblocks_y, get_block use ice_constants, only: c2 @@ -4809,8 +5407,9 @@ subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) !----------------------------------------------------------------------- integer (int_kind) :: & - i,j,iblk, &! dummy loop indices - numBlocks, &! number of local blocks + i,j,n,iblk,ii,jj, &! dummy loop indices + ilo,ihi,jlo,jhi, &! active block indices + numBlocks, &! number of local blocks blockID, &! block location ibc ! ghost cell column or row @@ -4831,32 +5430,40 @@ subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) do iblk = 1, numBlocks call ice_distributionGetBlockID(dist, iblk, blockID) this_block = get_block(blockID, blockID) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi if (this_block%iblock == 1) then ! west edge if (trim(ew_bndy_type) /= 'cyclic') then + do n = 1, nghost + ii = ilo - n ! gridcell to extrapolate to do j = 1, ny_block - ARRAY(1,j,iblk) = c2*ARRAY(2,j,iblk) - ARRAY(3,j,iblk) + ARRAY(ii,j,iblk) = c2*ARRAY(ii+1,j,iblk) - ARRAY(ii+2,j,iblk) + enddo enddo endif endif if (this_block%iblock == nblocks_x) then ! east edge if (trim(ew_bndy_type) /= 'cyclic') then - ! locate ghost cell column (avoid padding) - ibc = nx_block - do i = nx_block, nghost + 1, -1 - if (this_block%i_glob(i) == 0) ibc = ibc - 1 - enddo + do n = 1, nghost + ii = ihi + n ! gridcell to extrapolate to do j = 1, ny_block - ARRAY(ibc,j,iblk) = c2*ARRAY(ibc-1,j,iblk) - ARRAY(ibc-2,j,iblk) + ARRAY(ii,j,iblk) = c2*ARRAY(ii-1,j,iblk) - ARRAY(ii-2,j,iblk) + enddo enddo endif endif if (this_block%jblock == 1) then ! south edge if (trim(ns_bndy_type) /= 'cyclic') then + do n = 1, nghost + jj = jlo - n ! gridcell to extrapolate to do i = 1, nx_block - ARRAY(i,1,iblk) = c2*ARRAY(i,2,iblk) - ARRAY(i,3,iblk) + ARRAY(i,jj,iblk) = c2*ARRAY(i,jj+1,iblk) - ARRAY(i,jj+2,iblk) + enddo enddo endif endif @@ -4865,13 +5472,11 @@ subroutine ice_HaloExtrapolate2DR8(ARRAY,dist,ew_bndy_type,ns_bndy_type) if (trim(ns_bndy_type) /= 'cyclic' .and. & trim(ns_bndy_type) /= 'tripole' .and. & trim(ns_bndy_type) /= 'tripoleT' ) then - ! locate ghost cell column (avoid padding) - ibc = ny_block - do j = ny_block, nghost + 1, -1 - if (this_block%j_glob(j) == 0) ibc = ibc - 1 - enddo + do n = 1, nghost + jj = jhi + n ! gridcell to extrapolate to do i = 1, nx_block - ARRAY(i,ibc,iblk) = c2*ARRAY(i,ibc-1,iblk) - ARRAY(i,ibc-2,iblk) + ARRAY(i,jj,iblk) = c2*ARRAY(i,jj-1,iblk) - ARRAY(i,jj-2,iblk) + enddo enddo endif endif diff --git a/cicecore/cicedyn/infrastructure/comm/serial/ice_gather_scatter.F90 b/cicecore/cicedyn/infrastructure/comm/serial/ice_gather_scatter.F90 index 5f4938281..269183b97 100644 --- a/cicecore/cicedyn/infrastructure/comm/serial/ice_gather_scatter.F90 +++ b/cicecore/cicedyn/infrastructure/comm/serial/ice_gather_scatter.F90 @@ -925,80 +925,40 @@ subroutine scatter_global_dbl(ARRAY, ARRAY_G, src_task, dst_dist, & this_block = get_block(n,n) dst_block = dst_dist%blockLocalID(n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - ARRAY(i,j-yoffset2,dst_block) & - = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + ARRAY(i,j-yoffset2,dst_block) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + ARRAY(i,j,dst_block) = ARRAY_G(isrc,jsrc) + end do + endif + end do endif ! dst block not land end do ! block loop !----------------------------------------------------------------- - ! Ensure unused ghost cell values are 0 + ! Ensure/reset ghost cell values are 0 for noupdate !----------------------------------------------------------------- if (field_loc == field_loc_noupdate) then @@ -1173,80 +1133,40 @@ subroutine scatter_global_real(ARRAY, ARRAY_G, src_task, dst_dist, & this_block = get_block(n,n) dst_block = dst_dist%blockLocalID(n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - ARRAY(i,j-yoffset2,dst_block) & - = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + ARRAY(i,j-yoffset2,dst_block) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + ARRAY(i,j,dst_block) = ARRAY_G(isrc,jsrc) + end do + endif + end do endif ! dst block not land end do ! block loop !----------------------------------------------------------------- - ! Ensure unused ghost cell values are 0 + ! Ensure/reset ghost cell values are 0 for noupdate !----------------------------------------------------------------- if (field_loc == field_loc_noupdate) then @@ -1421,80 +1341,40 @@ subroutine scatter_global_int(ARRAY, ARRAY_G, src_task, dst_dist, & this_block = get_block(n,n) dst_block = dst_dist%blockLocalID(n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - ARRAY(i,j-yoffset2,dst_block) & - = isign * ARRAY_G(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + ARRAY(i,j-yoffset2,dst_block) & + = isign * ARRAY_G(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + ARRAY(i,j,dst_block) = ARRAY_G(isrc,jsrc) + end do + endif + end do endif ! dst block not land end do ! block loop !----------------------------------------------------------------- - ! Ensure unused ghost cell values are 0 + ! Ensure/reset ghost cell values are 0 for noupdate !----------------------------------------------------------------- if (field_loc == field_loc_noupdate) then @@ -1570,8 +1450,8 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) integer (int_kind) :: & i,j,n, &! dummy loop indices + isrc, jsrc, &! source addresses iblk, jblk, &! source addresses - iglb, jglb, &! global indices dst_block ! local block index in dest distribution type (block) :: & @@ -1618,13 +1498,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! southwest corner iblk = i jblk = j - iglb = this_block%i_glob(this_block%ilo)+i-1 - jglb = j - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ilo)+i-1 + jsrc = j + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) ! southeast corner iblk = this_block%ihi+i - iglb = this_block%i_glob(this_block%ihi)+nghost+i - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ihi)+nghost+i + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -1639,13 +1519,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northwest corner iblk = i jblk = this_block%jhi+j - iglb = this_block%i_glob(this_block%ilo)+i-1 - jglb = ny_global+nghost+j - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ilo)+i-1 + jsrc = ny_global+nghost+j + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) ! northeast corner iblk = this_block%ihi+i - iglb = this_block%i_glob(this_block%ihi)+nghost+i - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = this_block%i_glob(this_block%ihi)+nghost+i + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -1661,13 +1541,13 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northwest corner iblk = i jblk = this_block%jhi+j - iglb = i - jglb = this_block%j_glob(this_block%jhi)+nghost+j - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = i + jsrc = this_block%j_glob(this_block%jhi)+nghost+j + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) ! southwest corner jblk = j - jglb = this_block%j_glob(this_block%jlo)+j-1 - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + jsrc = this_block%j_glob(this_block%jlo)+j-1 + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) enddo enddo endif @@ -1684,17 +1564,16 @@ subroutine scatter_global_ext(ARRAY, ARRAY_G, src_task, dst_dist) ! northeast corner iblk = this_block%ihi+i jblk = this_block%jhi+j - iglb = nx_global+nghost+i - jglb = this_block%j_glob(this_block%jhi)+nghost+j - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + isrc = nx_global+nghost+i + jsrc = this_block%j_glob(this_block%jhi)+nghost+j + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) ! southeast corner jblk = j - jglb = this_block%j_glob(this_block%jlo)+j-1 - ARRAY(iblk,jblk,dst_block) = ARRAY_G(iglb,jglb) + jsrc = this_block%j_glob(this_block%jlo)+j-1 + ARRAY(iblk,jblk,dst_block) = ARRAY_G(isrc,jsrc) enddo enddo endif - endif ! dst block not land end do ! block loop @@ -1775,75 +1654,35 @@ subroutine scatter_global_stress(ARRAY, ARRAY_G1, ARRAY_G2, & this_block = get_block(n,n) dst_block = dst_dist%blockLocalID(n) - !*** if this is an interior block, then there is no - !*** padding or update checking required - - if (this_block%iblock > 1 .and. & - this_block%iblock < nblocks_x .and. & - this_block%jblock > 1 .and. & - this_block%jblock < nblocks_y) then - - do j=1,ny_block - do i=1,nx_block - ARRAY(i,j,dst_block) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - end do - end do - - !*** if this is an edge block but not a northern edge - !*** we only need to check for closed boundaries and - !*** padding (global index = 0) - - else if (this_block%jblock /= nblocks_y) then - - do j=1,ny_block - if (this_block%j_glob(j) /= 0) then - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - endif - end do - - !*** if this is a northern edge block, we need to check - !*** for and properly deal with tripole boundaries - - else - - do j=1,ny_block - if (this_block%j_glob(j) > 0) then ! normal boundary - - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - ARRAY(i,j,dst_block) = ARRAY_G1(this_block%i_glob(i),& - this_block%j_glob(j)) - endif - end do - - else if (this_block%j_glob(j) < 0) then ! tripole - - ! for yoffset=0 or 1, yoffset2=0,0 - ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid - do yoffset2=0,max(yoffset,0)-yoffset - jsrc = ny_global + yoffset + yoffset2 + & - (this_block%j_glob(j) + ny_global) - do i=1,nx_block - if (this_block%i_glob(i) /= 0) then - isrc = nx_global + xoffset - this_block%i_glob(i) - if (isrc < 1) isrc = isrc + nx_global - if (isrc > nx_global) isrc = isrc - nx_global - ARRAY(i,j-yoffset2,dst_block) & - = isign * ARRAY_G2(isrc,jsrc) - endif - end do - end do - - endif - end do - - endif + do j=1,ny_block + if (this_block%jblock == nblocks_y .and. this_block%j_glob(j) < 0) then + ! tripole is top block with j_glob < 0 + ! for yoffset=0 or 1, yoffset2=0,0 + ! for yoffset=-1, yoffset2=0,1, for u-rows on T-fold grid + do yoffset2=0,max(yoffset,0)-yoffset + jsrc = ny_global + yoffset + yoffset2 + & + (this_block%j_glob(j) + ny_global) + do i=1,nx_block + if (this_block%i_glob(i) /= 0) then + isrc = nx_global + xoffset - this_block%i_glob(i) + if (isrc < 1) isrc = isrc + nx_global + if (isrc > nx_global) isrc = isrc - nx_global + ARRAY(i,j-yoffset2,dst_block) & + = isign * ARRAY_G2(isrc,jsrc) + endif + end do + end do + else + ! normal block + do i=1,nx_block + isrc = this_block%i_glob(i) + jsrc = this_block%j_glob(j) + if (isrc >=1 .and. isrc <= nx_global .and. & + jsrc >=1 .and. jsrc <= ny_global) & + ARRAY(i,j,dst_block) = ARRAY_G1(isrc,jsrc) + end do + endif + end do endif ! dst block not land end do ! block loop diff --git a/cicecore/cicedyn/infrastructure/ice_blocks.F90 b/cicecore/cicedyn/infrastructure/ice_blocks.F90 index ccaf23999..de87c2a31 100644 --- a/cicecore/cicedyn/infrastructure/ice_blocks.F90 +++ b/cicecore/cicedyn/infrastructure/ice_blocks.F90 @@ -31,7 +31,13 @@ module ice_blocks tripoleTFlag ! tripole boundary is a T-fold integer (int_kind), dimension(:), pointer :: & - i_glob, j_glob ! global domain location for each point + i_glob, j_glob ! global domain location for each point. + ! valid values between 1:nx_global, 1:ny_global. + ! outside that range may occur in the halo with + ! open or closed bcs or on the tripole. + ! by definition, tripole is only on the north + ! boundary and in that case, the j_glob values + ! will be valid j_glob values with minus sign. end type public :: create_blocks ,& @@ -206,7 +212,7 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, & all_blocks_ij(iblock,jblock) = n do j=1,ny_block - j_global(j,n) = js - nghost + j - 1 + j_global(j,n) = js - nghost + j - 1 ! simple lower to upper counting !*** southern ghost cells @@ -215,13 +221,13 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, & case ('cyclic') j_global(j,n) = j_global(j,n) + ny_global case ('open') - j_global(j,n) = nghost - j + 1 + ! lower to upper case ('closed') - j_global(j,n) = 0 + ! lower to upper case ('tripole') - j_global(j,n) = nghost - j + 1 ! open + ! lower to upper case ('tripoleT') - j_global(j,n) = -j_global(j,n) + 1 ! open + ! lower to upper case default call abort_ice(subname//' ERROR: unknown n-s bndy type') end select @@ -239,13 +245,13 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, & case ('cyclic') j_global(j,n) = j_global(j,n) - ny_global case ('open') - j_global(j,n) = 2*ny_global - j_global(j,n) + 1 + ! lower to upper case ('closed') - j_global(j,n) = 0 + ! lower to upper case ('tripole') - j_global(j,n) = -j_global(j,n) + j_global(j,n) = -j_global(j,n) ! negative case ('tripoleT') - j_global(j,n) = -j_global(j,n) + j_global(j,n) = -j_global(j,n) ! negative case default call abort_ice(subname//' ERROR: unknown n-s bndy type') end select @@ -262,7 +268,7 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, & all_blocks(n)%j_glob => j_global(:,n) do i=1,nx_block - i_global(i,n) = is - nghost + i - 1 + i_global(i,n) = is - nghost + i - 1 ! left to right counting !*** western ghost cells @@ -271,9 +277,9 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, & case ('cyclic') i_global(i,n) = i_global(i,n) + nx_global case ('open') - i_global(i,n) = nghost - i + 1 + ! left to right case ('closed') - i_global(i,n) = 0 + ! left to right case default call abort_ice(subname//' ERROR: unknown e-w bndy type') end select @@ -291,9 +297,9 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, & case ('cyclic') i_global(i,n) = i_global(i,n) - nx_global case ('open') - i_global(i,n) = 2*nx_global - i_global(i,n) + 1 + ! left to right case ('closed') - i_global(i,n) = 0 + ! left to right case default call abort_ice(subname//' ERROR: unknown e-w bndy type') end select diff --git a/cicecore/cicedyn/infrastructure/ice_domain.F90 b/cicecore/cicedyn/infrastructure/ice_domain.F90 index 86d6a1939..e538d6b55 100644 --- a/cicecore/cicedyn/infrastructure/ice_domain.F90 +++ b/cicecore/cicedyn/infrastructure/ice_domain.F90 @@ -259,6 +259,12 @@ subroutine init_domain_blocks call abort_ice(subname//' ERROR: Not enough ghost cells allocated', file=__FILE__, line=__LINE__) endif + if (ns_boundary_type == 'closed') then + call abort_ice(subname//' ERROR: ns_boundary_type = '//trim(ns_boundary_type)//' not supported', file=__FILE__, line=__LINE__) + endif + if (ew_boundary_type == 'closed') then + call abort_ice(subname//' ERROR: ew_boundary_type = '//trim(ew_boundary_type)//' not supported', file=__FILE__, line=__LINE__) + endif !---------------------------------------------------------------------- ! ! compute block decomposition and details @@ -384,87 +390,6 @@ subroutine init_domain_distribution(KMTG,ULATG,grid_ice) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - if (trim(ns_boundary_type) == 'closed') then - call abort_ice(subname//' ERROR: ns_boundary_type = closed not supported', file=__FILE__, line=__LINE__) - allocate(nocn(nblocks_tot)) - nocn = 0 - do n=1,nblocks_tot - this_block = get_block(n,n) - if (this_block%jblock == nblocks_y) then ! north edge - do j = this_block%jhi-1, this_block%jhi - if (this_block%j_glob(j) > 0) then - do i = 1, nx_block - if (this_block%i_glob(i) > 0) then - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) - if (KMTG(ig,jg) > puny) nocn(n) = nocn(n) + 1 - endif - enddo - endif - enddo - endif - if (this_block%jblock == 1) then ! south edge - do j = this_block%jlo, this_block%jlo+1 - if (this_block%j_glob(j) > 0) then - do i = 1, nx_block - if (this_block%i_glob(i) > 0) then - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) - if (KMTG(ig,jg) > puny) nocn(n) = nocn(n) + 1 - endif - enddo - endif - enddo - endif - if (nocn(n) > 0) then - write(nu_diag,*) subname,'ns closed, Not enough land cells along ns edge' - call abort_ice(subname//' ERROR: Not enough land cells along ns edge for ns closed', & - file=__FILE__, line=__LINE__) - endif - enddo - deallocate(nocn) - endif - if (trim(ew_boundary_type) == 'closed') then - call abort_ice(subname//' ERROR: ew_boundary_type = closed not supported', file=__FILE__, line=__LINE__) - allocate(nocn(nblocks_tot)) - nocn = 0 - do n=1,nblocks_tot - this_block = get_block(n,n) - if (this_block%iblock == nblocks_x) then ! east edge - do j = 1, ny_block - if (this_block%j_glob(j) > 0) then - do i = this_block%ihi-1, this_block%ihi - if (this_block%i_glob(i) > 0) then - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) - if (KMTG(ig,jg) > puny) nocn(n) = nocn(n) + 1 - endif - enddo - endif - enddo - endif - if (this_block%iblock == 1) then ! west edge - do j = 1, ny_block - if (this_block%j_glob(j) > 0) then - do i = this_block%ilo, this_block%ilo+1 - if (this_block%i_glob(i) > 0) then - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) - if (KMTG(ig,jg) > puny) nocn(n) = nocn(n) + 1 - endif - enddo - endif - enddo - endif - if (nocn(n) > 0) then - write(nu_diag,*) subname,'ew closed, Not enough land cells along ew edge' - call abort_ice(subname//' ERROR: Not enough land cells along ew edge for ew closed', & - file=__FILE__, line=__LINE__) - endif - enddo - deallocate(nocn) - endif - !---------------------------------------------------------------------- ! ! estimate the amount of work per processor using the topography @@ -519,11 +444,11 @@ subroutine init_domain_distribution(KMTG,ULATG,grid_ice) do n=1,nblocks_tot this_block = get_block(n,n) do j=this_block%jlo,this_block%jhi - if (this_block%j_glob(j) > 0) then + jg = this_block%j_glob(j) + if (jg > 0) then do i=this_block%ilo,this_block%ihi - if (this_block%i_glob(i) > 0) then - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) + ig = this_block%i_glob(i) + if (ig > 0) then ! if (KMTG(ig,jg) > puny) & ! nocn(n) = max(nocn(n),nint(wght(ig,jg)+1.0_dbl_kind)) if (KMTG(ig,jg) > puny) then @@ -544,11 +469,11 @@ subroutine init_domain_distribution(KMTG,ULATG,grid_ice) do n=1,nblocks_tot this_block = get_block(n,n) do j=this_block%jlo,this_block%jhi - if (this_block%j_glob(j) > 0) then + jg = this_block%j_glob(j) + if (jg > 0) then do i=this_block%ilo,this_block%ihi - if (this_block%i_glob(i) > 0) then - ig = this_block%i_glob(i) - jg = this_block%j_glob(j) + ig = this_block%i_glob(i) + if (ig > 0) then if (grid_ice == 'C' .or. grid_ice == 'CD') then ! Have to be careful about block elimination with C/CD ! Use a bigger stencil diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index 4ae5448f5..e2ae52f6f 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -37,7 +37,7 @@ module ice_grid close_boundaries use ice_fileunits, only: nu_diag, nu_grid, nu_kmt, & get_fileunit, release_fileunit, flush_fileunit - use ice_gather_scatter, only: gather_global, scatter_global + use ice_gather_scatter, only: gather_global, scatter_global, gather_global_ext use ice_read_write, only: ice_read, ice_read_nc, ice_read_global, & ice_read_global_nc, ice_open, ice_open_nc, ice_close_nc, ice_check_nc use ice_timers, only: timer_bound, ice_timer_start, ice_timer_stop @@ -253,7 +253,62 @@ subroutine alloc_grid stat=ierr) if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory1', file=__FILE__, line=__LINE__) + dxT = c0 + dyT = c0 + dxU = c0 + dyU = c0 + dxN = c0 + dyN = c0 + dxE = c0 + dyE = c0 + HTE = c0 + HTN = c0 + tarea = c0 + uarea = c0 + narea = c0 + earea = c0 + tarear = c0 + uarear = c0 + narear = c0 + earear = c0 + tarean = c0 + tareas = c0 + ULON = c0 + ULAT = c0 + TLON = c0 + TLAT = c0 + NLON = c0 + NLAT = c0 + ELON = c0 + ELAT = c0 + ANGLE = c0 + ANGLET = c0 + bathymetry = c0 ocn_gridcell_frac(:,:,:) = -c1 ! special value to start, will be ignored unless set elsewhere + hm = c0 + bm = c0 + uvm = c0 + npm = c0 + epm = c0 + kmt = c0 + tmask = .false. + umask = .false. + umaskCD = .false. + nmask = .false. + emask = .false. + opmask = .false. + lmask_n = .false. + lmask_s = .false. + rndex_global = c0 + lont_bounds = c0 + latt_bounds = c0 + lonu_bounds = c0 + latu_bounds = c0 + lonn_bounds = c0 + latn_bounds = c0 + lone_bounds = c0 + late_bounds = c0 + if (save_ghte_ghtn) then if (my_task == master_task) then @@ -268,6 +323,8 @@ subroutine alloc_grid stat=ierr) endif if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory3', file=__FILE__, line=__LINE__) + G_HTE = c0 + G_HTN = c0 endif end subroutine alloc_grid @@ -595,14 +652,8 @@ subroutine init_grid2 !----------------------------------------------------------------- if (trim(grid_format) /= 'mom_nc') then - !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + !$OMP PARALLEL DO PRIVATE(iblk,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 = 1,ny_block do i = 1,nx_block tarea(i,j,iblk) = dxT(i,j,iblk)*dyT(i,j,iblk) @@ -615,13 +666,8 @@ subroutine init_grid2 !$OMP END PARALLEL DO endif - !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + !$OMP PARALLEL DO PRIVATE(iblk,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 = 1,ny_block do i = 1,nx_block @@ -959,6 +1005,8 @@ subroutine popgrid call ice_read_global(nu_grid,7,work_g1,'rda8',.true.) ! ANGLE call scatter_global(ANGLE, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_angle) + call ice_HaloExtrapolate(ANGLE, distrb_info, & + ew_boundary_type, ns_boundary_type) !----------------------------------------------------------------- ! cell dimensions @@ -1061,6 +1109,8 @@ subroutine popgrid_nc call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ANGLE call scatter_global(ANGLE, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_angle) + call ice_HaloExtrapolate(ANGLE, distrb_info, & + ew_boundary_type, ns_boundary_type) ! fix ANGLE: roundoff error due to single precision where (ANGLE > pi) ANGLE = pi where (ANGLE < -pi) ANGLE = -pi @@ -1081,16 +1131,22 @@ subroutine popgrid_nc call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) call scatter_global(ANGLET, work_g1, master_task, distrb_info, & field_loc_center, field_type_angle) + call ice_HaloExtrapolate(ANGLET, distrb_info, & + ew_boundary_type, ns_boundary_type) where (ANGLET > pi) ANGLET = pi where (ANGLET < -pi) ANGLET = -pi fieldname="tlon" call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) call scatter_global(TLON, work_g1, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(TLON, distrb_info, & + ew_boundary_type, ns_boundary_type) fieldname="tlat" call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) call scatter_global(TLAT, work_g1, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(TLAT, distrb_info, & + ew_boundary_type, ns_boundary_type) endif !----------------------------------------------------------------- ! cell dimensions @@ -1487,9 +1543,13 @@ subroutine mom_grid call mom_grid_rotation_angle(G_ULON, G_ULAT, G_TLON(1:nx_global,1:ny_global), work_g1) ! anglet call scatter_global(ANGLET, work_g1, master_task, distrb_info, & field_loc_center, field_type_angle) + call ice_HaloExtrapolate(ANGLET, distrb_info, & + ew_boundary_type, ns_boundary_type) call mom_grid_rotation_angle(G_TLON, G_TLAT, G_ULON(2:nx_global+1,2:ny_global+1), work_g1) ! angle call scatter_global(ANGLE, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_angle) + call ice_HaloExtrapolate(ANGLE, distrb_info, & + ew_boundary_type, ns_boundary_type) deallocate(work_g1, G_ULAT, G_TLAT, G_TLON, G_ULON, stat=ierr) if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) @@ -1792,26 +1852,28 @@ subroutine mom_dx(work_mom) jm1 = jm1 + 2 ; jm2 = jm2 + 2 enddo endif - - if (save_ghte_ghtn) then - do j = 1, ny_global - do i = 1, nx_global - G_HTN(i+nghost,j+nghost) = G_dxN(i,j) - enddo - enddo - call global_ext_halo(G_HTN) - endif endif call scatter_global(dxT, G_dxT, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dxT, distrb_info, & + ew_boundary_type, ns_boundary_type) call scatter_global(HTN, G_dxN, master_task, distrb_info, & field_loc_Nface, field_type_scalar) + call ice_HaloExtrapolate(HTN, distrb_info, & + ew_boundary_type, ns_boundary_type) + if (save_ghte_ghtn) then + call gather_global_ext(G_HTN, HTN, master_task, distrb_info) + endif dxN(:,:,:) = HTN(:,:,:) call scatter_global(dxE, G_dxE, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dxE, distrb_info, & + ew_boundary_type, ns_boundary_type) call scatter_global(dxU, G_dxU, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) + call ice_HaloExtrapolate(dxU, distrb_info, & + ew_boundary_type, ns_boundary_type) deallocate(G_dxT, G_dxE, G_dxU, G_dxN, stat=ierr) if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) @@ -1890,26 +1952,28 @@ subroutine mom_dy(work_mom) im1 = im1 + 2 ; im2 = im2 + 2 enddo endif - - if (save_ghte_ghtn) then - do j = 1, ny_global - do i = 1, nx_global - G_HTE(i+nghost,j+nghost) = G_dyE(i,j) - enddo - enddo - call global_ext_halo(G_HTE) - endif endif call scatter_global(dyT, G_dyT, master_task, distrb_info, & - field_loc_center, field_type_scalar) + field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dyT, distrb_info, & + ew_boundary_type, ns_boundary_type) call scatter_global(dyN, G_dyN, master_task, distrb_info, & field_loc_Nface, field_type_scalar) + call ice_HaloExtrapolate(dyN, distrb_info, & + ew_boundary_type, ns_boundary_type) call scatter_global(HTE, G_dyE, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(HTE, distrb_info, & + ew_boundary_type, ns_boundary_type) + if (save_ghte_ghtn) then + call gather_global_ext(G_HTE, HTE, master_task, distrb_info) + endif dyE(:,:,:) = HTE(:,:,:) call scatter_global(dyU, G_dyU, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) + call ice_HaloExtrapolate(dyU, distrb_info, & + ew_boundary_type, ns_boundary_type) deallocate(G_dyT, G_dyN, G_dyE, G_dyU) if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) @@ -2189,6 +2253,8 @@ subroutine geosgrid_nc call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) ! ANGLE call scatter_global(ANGLE, work_g1, master_task, distrb_info, & field_loc_NEcorner, field_type_angle) + call ice_HaloExtrapolate(ANGLE, distrb_info, & + ew_boundary_type, ns_boundary_type) ! fix ANGLE: roundoff error due to single precision where (ANGLE > pi) ANGLE = pi where (ANGLE < -pi) ANGLE = -pi @@ -2209,16 +2275,22 @@ subroutine geosgrid_nc call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) call scatter_global(ANGLET, work_g1, master_task, distrb_info, & field_loc_center, field_type_angle) + call ice_HaloExtrapolate(ANGLET, distrb_info, & + ew_boundary_type, ns_boundary_type) where (ANGLET > pi) ANGLET = pi where (ANGLET < -pi) ANGLET = -pi fieldname="tlon" call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) call scatter_global(TLON, work_g1, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(TLON, distrb_info, & + ew_boundary_type, ns_boundary_type) fieldname="tlat" call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag) call scatter_global(TLAT, work_g1, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(TLAT, distrb_info, & + ew_boundary_type, ns_boundary_type) endif !----------------------------------------------------------------- ! cell dimensions @@ -2358,6 +2430,10 @@ subroutine rectgrid call grid_boxislands_kmt(work_g1) + elseif (trim(kmt_type) == 'none') then + + work_g1(:,:) = c1 ! initialize hm as ocean + elseif (trim(kmt_type) == 'channel') then do j = 3,ny_global-2 ! closed top and bottom @@ -2597,7 +2673,6 @@ subroutine rectgrid_scale_dxdy call ice_HaloExtrapolate(ULAT, distrb_info, & ew_boundary_type, ns_boundary_type) - deallocate(work_g1) end subroutine rectgrid_scale_dxdy @@ -2775,19 +2850,18 @@ subroutine primary_grid_lengths_HTN(work_g) work_g2(i,j) = p5*(work_g(i,j) + work_g(ip1,j)) ! dxU enddo enddo - if (save_ghte_ghtn) then - do j = 1, ny_global - do i = 1,nx_global - G_HTN(i+nghost,j+nghost) = work_g(i,j) - enddo - enddo - call global_ext_halo(G_HTN) - endif endif call scatter_global(HTN, work_g, master_task, distrb_info, & field_loc_Nface, field_type_scalar) + call ice_HaloExtrapolate(HTN, distrb_info, & + ew_boundary_type, ns_boundary_type) + if (save_ghte_ghtn) then + call gather_global_ext(G_HTN, HTN, master_task, distrb_info) + endif call scatter_global(dxU, work_g2, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) + call ice_HaloExtrapolate(dxU, distrb_info, & + ew_boundary_type, ns_boundary_type) ! dxT = average of 2 neighbor HTNs in j @@ -2804,6 +2878,8 @@ subroutine primary_grid_lengths_HTN(work_g) endif call scatter_global(dxT, work_g2, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dxT, distrb_info, & + ew_boundary_type, ns_boundary_type) ! dxN = HTN @@ -2831,6 +2907,8 @@ subroutine primary_grid_lengths_HTN(work_g) endif call scatter_global(dxE, work_g2, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dxE, distrb_info, & + ew_boundary_type, ns_boundary_type) deallocate(work_g2, stat=ierr) if (ierr/=0) call abort_ice(subname//' ERROR: Dealloc error', file=__FILE__, line=__LINE__) @@ -2886,19 +2964,18 @@ subroutine primary_grid_lengths_HTE(work_g) work_g2(i,ny_global) = c2*work_g(i,ny_global-1) - work_g(i,ny_global-2) ! dyU enddo endif - if (save_ghte_ghtn) then - do j = 1, ny_global - do i = 1, nx_global - G_HTE(i+nghost,j+nghost) = work_g(i,j) - enddo - enddo - call global_ext_halo(G_HTE) - endif endif call scatter_global(HTE, work_g, master_task, distrb_info, & field_loc_Eface, field_type_scalar) + call ice_HaloExtrapolate(HTE, distrb_info, & + ew_boundary_type, ns_boundary_type) + if (save_ghte_ghtn) then + call gather_global_ext(G_HTE, HTE, master_task, distrb_info) + endif call scatter_global(dyU, work_g2, master_task, distrb_info, & field_loc_NEcorner, field_type_scalar) + call ice_HaloExtrapolate(dyU, distrb_info, & + ew_boundary_type, ns_boundary_type) ! dyT = average of 2 neighbor HTE in i @@ -2914,6 +2991,8 @@ subroutine primary_grid_lengths_HTE(work_g) endif call scatter_global(dyT, work_g2, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dyT, distrb_info, & + ew_boundary_type, ns_boundary_type) ! dyN = average of 4 neighbor HTEs @@ -2939,6 +3018,8 @@ subroutine primary_grid_lengths_HTE(work_g) endif call scatter_global(dyN, work_g2, master_task, distrb_info, & field_loc_center, field_type_scalar) + call ice_HaloExtrapolate(dyN, distrb_info, & + ew_boundary_type, ns_boundary_type) ! dyE = HTE diff --git a/cicecore/cicedyn/infrastructure/ice_read_write.F90 b/cicecore/cicedyn/infrastructure/ice_read_write.F90 index a6489473c..62545c3f3 100644 --- a/cicecore/cicedyn/infrastructure/ice_read_write.F90 +++ b/cicecore/cicedyn/infrastructure/ice_read_write.F90 @@ -1146,18 +1146,23 @@ subroutine ice_read_nc_xy(fid, nrec, varname, work, diag, & integer (kind=int_kind) :: lnrec ! local value of nrec - lnrec = nrec + logical (kind=log_kind) :: lrestart_ext ! local value of restart_ext - nx = nx_global - ny = ny_global + lnrec = nrec work = c0 ! to satisfy intent(out) attribute + lrestart_ext = .false. if (present(restart_ext)) then - if (restart_ext) then - nx = nx_global + 2*nghost - ny = ny_global + 2*nghost - endif + lrestart_ext = restart_ext + endif + + if (lrestart_ext) then + nx = nx_global + 2*nghost + ny = ny_global + 2*nghost + else + nx = nx_global + ny = ny_global endif if (my_task == master_task) then @@ -1243,10 +1248,8 @@ subroutine ice_read_nc_xy(fid, nrec, varname, work, diag, & ! NOTE: Ghost cells are not updated unless field_loc is present. !------------------------------------------------------------------- - if (present(restart_ext)) then - if (restart_ext) then - call scatter_global_ext(work, work_g1, master_task, distrb_info) - endif + if (lrestart_ext) then + call scatter_global_ext(work, work_g1, master_task, distrb_info) else if (present(field_loc)) then call scatter_global(work, work_g1, master_task, distrb_info, & @@ -1336,16 +1339,21 @@ subroutine ice_read_nc_xyz(fid, nrec, varname, work, diag, & integer (kind=int_kind) :: lnrec ! local value of nrec - lnrec = nrec + logical (kind=log_kind) :: lrestart_ext ! local value of restart_ext - nx = nx_global - ny = ny_global + lnrec = nrec + lrestart_ext = .false. if (present(restart_ext)) then - if (restart_ext) then - nx = nx_global + 2*nghost - ny = ny_global + 2*nghost - endif + lrestart_ext = restart_ext + endif + + if (lrestart_ext) then + nx = nx_global + 2*nghost + ny = ny_global + 2*nghost + else + nx = nx_global + ny = ny_global endif if (my_task == master_task) then @@ -1432,13 +1440,11 @@ subroutine ice_read_nc_xyz(fid, nrec, varname, work, diag, & ! NOTE: Ghost cells are not updated unless field_loc is present. !------------------------------------------------------------------- - if (present(restart_ext)) then - if (restart_ext) then - do n=1,ncat - call scatter_global_ext(work(:,:,n,:), work_g1(:,:,n), & - master_task, distrb_info) - enddo - endif + if (lrestart_ext) then + do n=1,ncat + call scatter_global_ext(work(:,:,n,:), work_g1(:,:,n), & + master_task, distrb_info) + enddo else if (present(field_loc)) then do n=1,ncat @@ -1532,20 +1538,25 @@ subroutine ice_read_nc_xyf(fid, nrec, varname, work, diag, & integer (kind=int_kind) :: lnrec ! local value of nrec + logical (kind=log_kind) :: lrestart_ext ! local value of restart_ext + character(len=*), parameter :: subname = '(ice_read_nc_xyf)' #ifdef USE_NETCDF lnrec = nrec - nx = nx_global - ny = ny_global - + lrestart_ext = .false. if (present(restart_ext)) then - if (restart_ext) then - nx = nx_global + 2*nghost - ny = ny_global + 2*nghost - endif + lrestart_ext = restart_ext + endif + + if (lrestart_ext) then + nx = nx_global + 2*nghost + ny = ny_global + 2*nghost + else + nx = nx_global + ny = ny_global endif if (my_task == master_task) then @@ -1632,13 +1643,11 @@ subroutine ice_read_nc_xyf(fid, nrec, varname, work, diag, & ! NOTE: Ghost cells are not updated unless field_loc is present. !------------------------------------------------------------------- - if (present(restart_ext)) then - if (restart_ext) then - do n = 1, nfreq - call scatter_global_ext(work(:,:,n,1,:), work_g1(:,:,n), & - master_task, distrb_info) - enddo - endif + if (lrestart_ext) then + do n = 1, nfreq + call scatter_global_ext(work(:,:,n,1,:), work_g1(:,:,n), & + master_task, distrb_info) + enddo else if (present(field_loc)) then do n = 1, nfreq @@ -2214,14 +2223,19 @@ subroutine ice_write_nc_xy(fid, nrec, varid, work, diag, & integer (kind=int_kind) :: nx, ny - nx = nx_global - ny = ny_global + logical (kind=log_kind) :: lrestart_ext ! local value of restart_ext + lrestart_ext = .false. if (present(restart_ext)) then - if (restart_ext) then - nx = nx_global + 2*nghost - ny = ny_global + 2*nghost - endif + lrestart_ext = restart_ext + endif + + if (lrestart_ext) then + nx = nx_global + 2*nghost + ny = ny_global + 2*nghost + else + nx = nx_global + ny = ny_global endif if (present(varname)) then @@ -2236,10 +2250,8 @@ subroutine ice_write_nc_xy(fid, nrec, varid, work, diag, & allocate(work_g1(1,1)) ! to save memory endif - if (present(restart_ext)) then - if (restart_ext) then - call gather_global_ext(work_g1, work, master_task, distrb_info, spc_val=c0) - endif + if (lrestart_ext) then + call gather_global_ext(work_g1, work, master_task, distrb_info, spc_val=c0) else call gather_global(work_g1, work, master_task, distrb_info, spc_val=c0) endif @@ -2338,14 +2350,19 @@ subroutine ice_write_nc_xyz(fid, nrec, varid, work, diag, & integer (kind=int_kind) :: nx, ny - nx = nx_global - ny = ny_global + logical (kind=log_kind) :: lrestart_ext ! local value of restart_ext + lrestart_ext = .false. if (present(restart_ext)) then - if (restart_ext) then - nx = nx_global + 2*nghost - ny = ny_global + 2*nghost - endif + lrestart_ext = restart_ext + endif + + if (lrestart_ext) then + nx = nx_global + 2*nghost + ny = ny_global + 2*nghost + else + nx = nx_global + ny = ny_global endif if (my_task == master_task) then @@ -2354,13 +2371,11 @@ subroutine ice_write_nc_xyz(fid, nrec, varid, work, diag, & allocate(work_g1(1,1,ncat)) ! to save memory endif - if (present(restart_ext)) then - if (restart_ext) then - do n=1,ncat - call gather_global_ext(work_g1(:,:,n), work(:,:,n,:), & - master_task, distrb_info, spc_val=c0) - enddo - endif + if (lrestart_ext) then + do n=1,ncat + call gather_global_ext(work_g1(:,:,n), work(:,:,n,:), & + master_task, distrb_info, spc_val=c0) + enddo else do n=1,ncat call gather_global(work_g1(:,:,n), work(:,:,n,:), & @@ -2664,14 +2679,19 @@ subroutine ice_read_nc_uv(fid, nrec, nzlev, varname, work, diag, & integer (kind=int_kind) :: nx, ny - nx = nx_global - ny = ny_global + logical (kind=log_kind) :: lrestart_ext ! local value of restart_ext + lrestart_ext = .false. if (present(restart_ext)) then - if (restart_ext) then - nx = nx_global + 2*nghost - ny = ny_global + 2*nghost - endif + lrestart_ext = restart_ext + endif + + if (lrestart_ext) then + nx = nx_global + 2*nghost + ny = ny_global + 2*nghost + else + nx = nx_global + ny = ny_global endif if (my_task == master_task) then @@ -2717,10 +2737,8 @@ subroutine ice_read_nc_uv(fid, nrec, nzlev, varname, work, diag, & ! NOTE: Ghost cells are not updated unless field_loc is present. !------------------------------------------------------------------- - if (present(restart_ext)) then - if (restart_ext) then - call scatter_global_ext(work, work_g1, master_task, distrb_info) - endif + if (lrestart_ext) then + call scatter_global_ext(work, work_g1, master_task, distrb_info) else if (present(field_loc)) then call scatter_global(work, work_g1, master_task, distrb_info, & diff --git a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_restart.F90 b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_restart.F90 index db6aa80bf..1c25e3f30 100644 --- a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_restart.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_restart.F90 @@ -10,6 +10,7 @@ module ice_restart use ice_broadcast + use ice_constants, only: c0 use ice_communicate, only: my_task, master_task use ice_kinds_mod #ifdef USE_NETCDF @@ -748,39 +749,25 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3, & character(len=*), parameter :: subname = '(read_restart_field)' + work (:,:,:,:) = c0 + work2(:,:,:) = c0 #ifdef USE_NETCDF if (present(field_loc)) then if (ndim3 == ncat) then - if (restart_ext) then - call ice_read_nc(ncid,1,vname,work,diag, & - field_loc=field_loc,field_type=field_type,restart_ext=restart_ext) - else - call ice_read_nc(ncid,1,vname,work,diag,field_loc,field_type) - endif + call ice_read_nc(ncid,1,vname,work,diag, & + field_loc=field_loc,field_type=field_type,restart_ext=restart_ext) elseif (ndim3 == 1) then - if (restart_ext) then - call ice_read_nc(ncid,1,vname,work2,diag, & - field_loc=field_loc,field_type=field_type,restart_ext=restart_ext) - else - call ice_read_nc(ncid,1,vname,work2,diag,field_loc,field_type) - endif + call ice_read_nc(ncid,1,vname,work2,diag, & + field_loc=field_loc,field_type=field_type,restart_ext=restart_ext) work(:,:,1,:) = work2(:,:,:) else write(nu_diag,*) 'ndim3 not supported ',ndim3 endif else if (ndim3 == ncat) then - if (restart_ext) then - call ice_read_nc(ncid, 1, vname, work, diag, restart_ext=restart_ext) - else - call ice_read_nc(ncid, 1, vname, work, diag) - endif + call ice_read_nc(ncid, 1, vname, work, diag, restart_ext=restart_ext) elseif (ndim3 == 1) then - if (restart_ext) then - call ice_read_nc(ncid, 1, vname, work2, diag, restart_ext=restart_ext) - else - call ice_read_nc(ncid, 1, vname, work2, diag) - endif + call ice_read_nc(ncid, 1, vname, work2, diag, restart_ext=restart_ext) work(:,:,1,:) = work2(:,:,:) else write(nu_diag,*) 'ndim3 not supported ',ndim3 @@ -841,18 +828,10 @@ subroutine write_restart_field(nu,nrec,work,atype,vname,ndim3,diag) call ice_check_nc(status, subname//' ERROR: inq varid '//trim(vname), file=__FILE__, line=__LINE__) endif if (ndim3 == ncat) then - if (restart_ext) then - call ice_write_nc(ncid, 1, varid, work, diag, restart_ext, varname=trim(vname)) - else - call ice_write_nc(ncid, 1, varid, work, diag, varname=trim(vname)) - endif + call ice_write_nc(ncid, 1, varid, work, diag, restart_ext, varname=trim(vname)) elseif (ndim3 == 1) then work2(:,:,:) = work(:,:,1,:) - if (restart_ext) then - call ice_write_nc(ncid, 1, varid, work2, diag, restart_ext, varname=trim(vname)) - else - call ice_write_nc(ncid, 1, varid, work2, diag, varname=trim(vname)) - endif + call ice_write_nc(ncid, 1, varid, work2, diag, restart_ext, varname=trim(vname)) else write(nu_diag,*) 'ndim3 not supported',ndim3 endif diff --git a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 index 11f0fb803..37cf4d985 100644 --- a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_restart.F90 @@ -756,6 +756,7 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3,diag, & subname// " ERROR: missing varndims "//trim(vname),file=__FILE__,line=__LINE__) call pio_seterrorhandling(File, PIO_INTERNAL_ERROR) + work (:,:,:,:) = c0 if (ndim3 == ncat .and. ndims == 3) then call pio_read_darray(File, vardesc, iodesc3d_ncat, work, status) #ifdef CESMCOUPLED diff --git a/cicecore/drivers/unittest/gridavgchk/gridavgchk.F90 b/cicecore/drivers/unittest/gridavgchk/gridavgchk.F90 index bd7ed3165..02e51fbba 100644 --- a/cicecore/drivers/unittest/gridavgchk/gridavgchk.F90 +++ b/cicecore/drivers/unittest/gridavgchk/gridavgchk.F90 @@ -16,7 +16,6 @@ program gridavgchk use CICE_InitMod use ice_kinds_mod, only: int_kind, dbl_kind use ice_blocks, only: block, get_block, nx_block, ny_block, nblocks_tot - use ice_boundary, only: ice_HaloUpdate use ice_constants, only: c0, c1, c2, p25, & field_loc_center, field_loc_NEcorner, & field_loc_Nface, field_loc_Eface, field_type_scalar diff --git a/cicecore/drivers/unittest/halochk/halochk.F90 b/cicecore/drivers/unittest/halochk/halochk.F90 index 29eaa8150..2fecaf31e 100644 --- a/cicecore/drivers/unittest/halochk/halochk.F90 +++ b/cicecore/drivers/unittest/halochk/halochk.F90 @@ -40,11 +40,12 @@ program halochk implicit none - integer(int_kind) :: nn, nl, nt, i, j, k1, k2, n, ib, ie, jb, je + integer(int_kind) :: nn, nl, nt, nf, i, j, k1, k2, n, ib, ie, jb, je integer(int_kind) :: iblock, itrip, ioffset, joffset integer(int_kind) :: blockID, numBlocks, jtrip type (block) :: this_block + ! fields sent to the haloupdate real(dbl_kind) , allocatable :: darrayi1(:,:,:) , darrayj1(:,:,:) real(dbl_kind) , allocatable :: darrayi2(:,:,:,:) , darrayj2(:,:,:,:) real(dbl_kind) , allocatable :: darrayi3(:,:,:,:,:), darrayj3(:,:,:,:,:) @@ -58,25 +59,27 @@ program halochk real(dbl_kind) , allocatable :: darrayi1str(:,:,:) , darrayj1str(:,:,:) real(dbl_kind) , allocatable :: darrayi10(:,:,:) , darrayj10(:,:,:) - real(dbl_kind), allocatable :: cidata_bas(:,:,:,:,:),cjdata_bas(:,:,:,:,:) - real(dbl_kind), allocatable :: cidata_nup(:,:,:,:,:),cjdata_nup(:,:,:,:,:) - real(dbl_kind), allocatable :: cidata_std(:,:,:,:,:),cjdata_std(:,:,:,:,:) + ! expected results + real(dbl_kind), allocatable :: cidata_bas(:,:,:,:,:),cjdata_bas(:,:,:,:,:) ! baseline integer(int_kind), parameter :: maxtests = 11 integer(int_kind), parameter :: maxtypes = 4 integer(int_kind), parameter :: maxlocs = 5 + integer(int_kind), parameter :: maxfills = 2 integer(int_kind), parameter :: nz1 = 3 integer(int_kind), parameter :: nz2 = 4 real(dbl_kind) :: aichk,ajchk,cichk,cjchk,rival,rjval,rsign - character(len=16) :: locs_name(maxlocs), types_name(maxtypes) + real(dbl_kind) :: fillexpected + character(len=16) :: locs_name(maxlocs), types_name(maxtypes), fill_name(maxfills) integer(int_kind) :: field_loc(maxlocs), field_type(maxtypes) + logical :: halofill integer(int_kind) :: npes, ierr, ntask, testcnt, tottest, tpcnt, tfcnt integer(int_kind) :: errorflag0, gflag, k1m, k2m, ptcntsum, failcntsum integer(int_kind), allocatable :: errorflag(:) integer(int_kind), allocatable :: ptcnt(:), failcnt(:) character(len=128), allocatable :: teststring(:) character(len=32) :: halofld - logical :: tripole_average, tripole_pole, spvalL1 + logical :: tripole_average, tripole_pole logical :: first_call = .true. real(dbl_kind) , parameter :: fillval = -88888.0_dbl_kind @@ -94,6 +97,7 @@ program halochk locs_name (:) = 'unknown' types_name(:) = 'unknown' + fill_name (:) = 'unknown' field_type(:) = field_type_unknown field_loc (:) = field_loc_unknown @@ -110,7 +114,7 @@ program halochk locs_name (1) = 'center' field_loc (1) = field_loc_center - locs_name (2) = 'NEcorner' + locs_name (2) = 'NEcorn' field_loc (2) = field_loc_NEcorner locs_name (3) = 'Nface' field_loc (3) = field_loc_Nface @@ -121,7 +125,10 @@ program halochk ! locs_name (6) = 'unknown' ! field_loc (6) = field_loc_unknown ! aborts in CICE, as expected - tottest = maxtests * maxlocs * maxtypes + fill_name (1) = 'fill' + fill_name (2) = 'nofill' + + tottest = maxtests * maxlocs * maxtypes * maxfills allocate(errorflag(tottest)) allocate(teststring(tottest)) allocate(ptcnt(tottest)) @@ -187,10 +194,6 @@ program halochk allocate(cidata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) allocate(cjdata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) - allocate(cidata_std(nx_block,ny_block,nz1,nz2,max_blocks)) - allocate(cjdata_std(nx_block,ny_block,nz1,nz2,max_blocks)) - allocate(cidata_nup(nx_block,ny_block,nz1,nz2,max_blocks)) - allocate(cjdata_nup(nx_block,ny_block,nz1,nz2,max_blocks)) darrayi1 = fillval darrayj1 = fillval @@ -218,14 +221,14 @@ program halochk darrayj10 = fillval cidata_bas = fillval cjdata_bas = fillval - cidata_std = fillval - cjdata_std = fillval - cidata_nup = fillval - cjdata_nup = fillval call ice_distributionGet(distrb_info, numLocalBlocks = numBlocks) !--- baseline data --- + ! set to the global index + ! i/j valid everywhere for "cyclic" + ! i/j valid for "open" with extrapolation on outer boundary + ! i/j zero on outer boundary for "closed" do iblock = 1,numBlocks call ice_distributionGetBlockID(distrb_info, iblock, blockID) @@ -244,102 +247,32 @@ program halochk enddo enddo - !--- setup nup (noupdate) solution, set halo/pad will fillval --- - - cidata_nup(:,:,:,:,:) = cidata_bas(:,:,:,:,:) - cjdata_nup(:,:,:,:,:) = cjdata_bas(:,:,:,:,:) - - do iblock = 1,numBlocks - call ice_distributionGetBlockID(distrb_info, iblock, blockID) - this_block = get_block(blockID, blockID) - ib = this_block%ilo - ie = this_block%ihi - jb = this_block%jlo - je = this_block%jhi - cidata_nup(1:ib-1 ,: ,:,:,iblock) = fillval - cjdata_nup(1:ib-1 ,: ,:,:,iblock) = fillval - cidata_nup(ie+1:nx_block,: ,:,:,iblock) = fillval - cjdata_nup(ie+1:nx_block,: ,:,:,iblock) = fillval - cidata_nup(: ,1:jb-1 ,:,:,iblock) = fillval - cjdata_nup(: ,1:jb-1 ,:,:,iblock) = fillval - cidata_nup(: ,je+1:ny_block,:,:,iblock) = fillval - cjdata_nup(: ,je+1:ny_block,:,:,iblock) = fillval - enddo - - !--- setup std solution for cyclic, closed, open, tripole solution --- - - cidata_std(:,:,:,:,:) = cidata_bas(:,:,:,:,:) - cjdata_std(:,:,:,:,:) = cjdata_bas(:,:,:,:,:) - - !--- halo off on east and west boundary --- - if (ew_boundary_type == 'closed' .or. & - ew_boundary_type == 'open' ) then - do iblock = 1,numBlocks - call ice_distributionGetBlockID(distrb_info, iblock, blockID) - this_block = get_block(blockID, blockID) - ib = this_block%ilo - ie = this_block%ihi - jb = this_block%jlo - je = this_block%jhi - if (this_block%i_glob(ib) == 1) then - cidata_std(1:ib-1 ,:,:,:,iblock) = dhalofillval - cjdata_std(1:ib-1 ,:,:,:,iblock) = dhalofillval - endif - if (this_block%i_glob(ie) == nx_global) then - cidata_std(ie+1:nx_block,:,:,:,iblock) = dhalofillval - cjdata_std(ie+1:nx_block,:,:,:,iblock) = dhalofillval - endif - enddo - endif - - !--- halo off on south boundary --- - if (ns_boundary_type == 'closed' .or. & - ns_boundary_type == 'open' .or. & - ns_boundary_type == 'tripole' .or. & - ns_boundary_type == 'tripoleT' ) then - do iblock = 1,numBlocks - call ice_distributionGetBlockID(distrb_info, iblock, blockID) - this_block = get_block(blockID, blockID) - ib = this_block%ilo - ie = this_block%ihi - jb = this_block%jlo - je = this_block%jhi - if (this_block%j_glob(jb) == 1) then - cidata_std(:,1:jb-1,:,:,iblock) = dhalofillval - cjdata_std(:,1:jb-1,:,:,iblock) = dhalofillval - endif - enddo - endif - - !--- halo off on north boundary, tripole handled later --- - if (ns_boundary_type == 'closed' .or. & - ns_boundary_type == 'open' .or. & - ns_boundary_type == 'tripole' .or. & - ns_boundary_type == 'tripoleT' ) then - do iblock = 1,numBlocks - call ice_distributionGetBlockID(distrb_info, iblock, blockID) - this_block = get_block(blockID, blockID) - ib = this_block%ilo - ie = this_block%ihi - jb = this_block%jlo - je = this_block%jhi - if (this_block%j_glob(je) == ny_global) then - cidata_std(:,je+1:ny_block,:,:,iblock) = dhalofillval - cjdata_std(:,je+1:ny_block,:,:,iblock) = dhalofillval - endif - enddo - endif - !--------------------------------------------------------------- testcnt = 0 do nn = 1, maxtests do nl = 1, maxlocs do nt = 1, maxtypes + do nf = 1, maxfills !--- setup test --- first_call = .true. testcnt = testcnt + 1 + if (nf == 1) then + halofill = .true. + fillexpected = dhalofillval + elseif (nf == 2) then + halofill = .false. + fillexpected = fillval + else + write(6,*) subname,' nf = ',nf + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) 'HALOCHK FAILED' + write(6,*) ' ' + endif + call abort_ice(subname//' invalid value of nf',file=__FILE__,line=__LINE__) + endif if (testcnt > tottest) then if (my_task == master_task) then write(6,*) ' ' @@ -388,35 +321,54 @@ program halochk darrayi10 = darrayi1 darrayj10 = darrayj1 - !--- halo update --- if (nn == 1) then k1m = 1 k2m = 1 halofld = '2DR8' - call ice_haloUpdate(darrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) - call ice_haloUpdate(darrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + if (halofill) then + call ice_haloUpdate(darrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + else + call ice_haloUpdate(darrayi1, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(darrayj1, halo_info, field_loc(nl), field_type(nt)) + endif elseif (nn == 2) then k1m = nz1 k2m = 1 halofld = '3DR8' - call ice_haloUpdate(darrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) - call ice_haloUpdate(darrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + if (halofill) then + call ice_haloUpdate(darrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + else + call ice_haloUpdate(darrayi2, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(darrayj2, halo_info, field_loc(nl), field_type(nt)) + endif elseif (nn == 3) then k1m = nz1 k2m = nz2 halofld = '4DR8' - call ice_haloUpdate(darrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) - call ice_haloUpdate(darrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + if (halofill) then + call ice_haloUpdate(darrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + else + call ice_haloUpdate(darrayi3, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(darrayj3, halo_info, field_loc(nl), field_type(nt)) + endif elseif (nn == 4) then k1m = 1 k2m = 1 halofld = '2DR4' rarrayi1 = real(darrayi1,kind=real_kind) rarrayj1 = real(darrayj1,kind=real_kind) - call ice_haloUpdate(rarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) - call ice_haloUpdate(rarrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + if (halofill) then + call ice_haloUpdate(rarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + else + call ice_haloUpdate(rarrayi1, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(rarrayj1, halo_info, field_loc(nl), field_type(nt)) + endif darrayi1 = real(rarrayi1,kind=dbl_kind) darrayj1 = real(rarrayj1,kind=dbl_kind) elseif (nn == 5) then @@ -425,8 +377,13 @@ program halochk halofld = '3DR4' rarrayi2 = real(darrayi2,kind=real_kind) rarrayj2 = real(darrayj2,kind=real_kind) - call ice_haloUpdate(rarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) - call ice_haloUpdate(rarrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + if (halofill) then + call ice_haloUpdate(rarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + else + call ice_haloUpdate(rarrayi2, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(rarrayj2, halo_info, field_loc(nl), field_type(nt)) + endif darrayi2 = real(rarrayi2,kind=dbl_kind) darrayj2 = real(rarrayj2,kind=dbl_kind) elseif (nn == 6) then @@ -435,8 +392,13 @@ program halochk halofld = '4DR4' rarrayi3 = real(darrayi3,kind=real_kind) rarrayj3 = real(darrayj3,kind=real_kind) - call ice_haloUpdate(rarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) - call ice_haloUpdate(rarrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + if (halofill) then + call ice_haloUpdate(rarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + else + call ice_haloUpdate(rarrayi3, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(rarrayj3, halo_info, field_loc(nl), field_type(nt)) + endif darrayi3 = real(rarrayi3,kind=dbl_kind) darrayj3 = real(rarrayj3,kind=dbl_kind) elseif (nn == 7) then @@ -445,8 +407,13 @@ program halochk halofld = '2DI4' iarrayi1 = nint(darrayi1) iarrayj1 = nint(darrayj1) - call ice_haloUpdate(iarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) - call ice_haloUpdate(iarrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + if (halofill) then + call ice_haloUpdate(iarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + else + call ice_haloUpdate(iarrayi1, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(iarrayj1, halo_info, field_loc(nl), field_type(nt)) + endif darrayi1 = real(iarrayi1,kind=dbl_kind) darrayj1 = real(iarrayj1,kind=dbl_kind) elseif (nn == 8) then @@ -455,8 +422,13 @@ program halochk halofld = '3DI4' iarrayi2 = nint(darrayi2) iarrayj2 = nint(darrayj2) - call ice_haloUpdate(iarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) - call ice_haloUpdate(iarrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + if (halofill) then + call ice_haloUpdate(iarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + else + call ice_haloUpdate(iarrayi2, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(iarrayj2, halo_info, field_loc(nl), field_type(nt)) + endif darrayi2 = real(iarrayi2,kind=dbl_kind) darrayj2 = real(iarrayj2,kind=dbl_kind) elseif (nn == 9) then @@ -465,20 +437,36 @@ program halochk halofld = '4DI4' iarrayi3 = nint(darrayi3) iarrayj3 = nint(darrayj3) - call ice_haloUpdate(iarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) - call ice_haloUpdate(iarrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + if (halofill) then + call ice_haloUpdate(iarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + else + call ice_haloUpdate(iarrayi3, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(iarrayj3, halo_info, field_loc(nl), field_type(nt)) + endif darrayi3 = real(iarrayi3,kind=dbl_kind) darrayj3 = real(iarrayj3,kind=dbl_kind) elseif (nn == 10) then k1m = 1 k2m = 1 halofld = '2DL1' - larrayi1 = .true. - where (darrayi1 == fillval) larrayi1 = .false. - larrayj1 = .false. - where (darrayj1 == fillval) larrayj1 = .true. - call ice_haloUpdate(larrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=0) - call ice_haloUpdate(larrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=1) + where (darrayi1 == fillval) + larrayi1 = .false. + elsewhere + larrayi1 = (mod(nint(darrayi1),2) == 1) + endwhere + where (darrayj1 == fillval) + larrayj1 = .true. + elsewhere + larrayj1 = (mod(nint(darrayj1),2) == 1) + endwhere + if (halofill) then + call ice_haloUpdate(larrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=0) + call ice_haloUpdate(larrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=1) + else + call ice_haloUpdate(larrayi1, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate(larrayj1, halo_info, field_loc(nl), field_type(nt)) + endif darrayi1 = c0 where (larrayi1) darrayi1 = c1 darrayj1 = c0 @@ -489,11 +477,16 @@ program halochk halofld = 'STRESS' darrayi1str = -darrayi1 ! flip sign for testing darrayj1str = -darrayj1 - call ice_haloUpdate_stress(darrayi1, darrayi1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) - call ice_haloUpdate_stress(darrayj1, darrayj1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + if (halofill) then + call ice_haloUpdate_stress(darrayi1, darrayi1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate_stress(darrayj1, darrayj1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + else + call ice_haloUpdate_stress(darrayi1, darrayi1str, halo_info, field_loc(nl), field_type(nt)) + call ice_haloUpdate_stress(darrayj1, darrayj1str, halo_info, field_loc(nl), field_type(nt)) + endif endif - write(teststring(testcnt),'(5a10)') trim(halofld),trim(locs_name(nl)),trim(types_name(nt)), & + write(teststring(testcnt),'(6a8)') trim(halofld),trim(locs_name(nl)),trim(types_name(nt)),trim(fill_name(nf)), & trim(ew_boundary_type),trim(ns_boundary_type) do iblock = 1,numBlocks @@ -504,15 +497,12 @@ program halochk jb = this_block%jlo je = this_block%jhi ! just check non-padded gridcells -! do j = 1,ny_block -! do i = 1,nx_block do j = jb-nghost, je+nghost do i = ib-nghost, ie+nghost do k1 = 1,k1m do k2 = 1,k2m tripole_average = .false. tripole_pole = .false. - spvalL1 = .false. if (index(halofld,'2D') > 0) then aichk = darrayi1(i,j,iblock) ajchk = darrayj1(i,j,iblock) @@ -534,14 +524,46 @@ program halochk call abort_ice(subname//' halofld not matched '//trim(halofld),file=__FILE__,line=__LINE__) endif + cichk = cidata_bas(i,j,k1,k2,iblock) + cjchk = cjdata_bas(i,j,k1,k2,iblock) + + ! halo special cases if (field_loc (nl) == field_loc_noupdate .or. & field_type(nt) == field_type_noupdate) then - cichk = cidata_nup(i,j,k1,k2,iblock) - cjchk = cjdata_nup(i,j,k1,k2,iblock) + if (i < ib .or. j < jb .or. i > ie .or. j > je) then + ! no halo update anywhere, doesn't even see fillvalue passed in + cichk = fillval + cjchk = fillval + endif + else - cichk = cidata_std(i,j,k1,k2,iblock) - cjchk = cjdata_std(i,j,k1,k2,iblock) + ! if ew_boundary_type is not cyclic we expect just fill values on outer boundary + if (ew_boundary_type /= 'cyclic' .and. & + ((this_block%i_glob(ib) == 1 .and. i < ib) .or. & ! west outer face + (this_block%i_glob(ie) == nx_global .and. i > ie))) then ! east outer face + cichk = fillexpected + cjchk = fillexpected + endif + + ! if ns_boundary_type is not cyclic we expect just fill values on outer boundary except + ! - tripole north edge will be haloed and is updated below, default to fill value for now + ! - tripole south edge will be set to the fillvalue or to haloupdate internal default (c0) + ! tripole basically assumes south edge is land or always ice free in CICE + if (ns_boundary_type /= 'cyclic' .and. & + ((this_block%j_glob(jb) == 1 .and. j < jb) .or. & ! south outer face + (this_block%j_glob(je) == ny_global .and. j > je))) then ! north outer face + ! ns_boundary_type is not cyclic and on outer boundary + if ((ns_boundary_type == 'tripole' .or. & + ns_boundary_type == 'tripoleT') .and. & + .not. halofill) then + cichk = c0 + cjchk = c0 + else + cichk = fillexpected + cjchk = fillexpected + endif + endif if (index(halofld,'STRESS') > 0) then ! only updates on tripole zipper for tripole grids @@ -560,11 +582,11 @@ program halochk (ns_boundary_type == 'tripoleT' .and. & (j >= je)))) then - ! flip sign for vector/angle - if (field_type(nt) == field_type_vector .or. field_type(nt) == field_type_angle ) then + ! flip sign for vector/angle except for logical halo updates + rsign = c1 + if ((field_type(nt) == field_type_vector .or. field_type(nt) == field_type_angle) .and. & + .not. (index(halofld,'L1') > 0)) then rsign = -c1 - else - rsign = c1 endif ! for tripole @@ -650,44 +672,40 @@ program halochk if (index(halofld,'STRESS') > 0) then ! only updates on tripole zipper for tripole grids, not tripoleT + ! note: L1 and STRESS never overlap so don't worry about L1 here if (tripole_pole) then ! flip sign due to sign of darrayi1str ! ends of tripole seam not averaged in CICE - cichk = -rsign * cidata_std(i,j,k1,k2,iblock) - cjchk = -rsign * cjdata_std(i,j,k1,k2,iblock) + cichk = -rsign * cidata_bas(i,j,k1,k2,iblock) + cjchk = -rsign * cjdata_bas(i,j,k1,k2,iblock) else cichk = -rsign * rival cjchk = -rsign * rjval endif - elseif (index(halofld,'L1') > 0 .and. j == je) then - ! force cichk and cjchk to match on tripole average index, calc not well defined - spvalL1 = .true. - cichk = aichk - cjchk = ajchk + elseif (tripole_pole) then ! ends of tripole seam not averaged in CICE - cichk = rsign * cidata_std(i,j,k1,k2,iblock) - cjchk = rsign * cjdata_std(i,j,k1,k2,iblock) + cichk = rsign * cidata_bas(i,j,k1,k2,iblock) + cjchk = rsign * cjdata_bas(i,j,k1,k2,iblock) + elseif (tripole_average) then - ! tripole average - cichk = p5 * (cidata_std(i,j,k1,k2,iblock) + rsign * rival) - cjchk = p5 * (cjdata_std(i,j,k1,k2,iblock) + rsign * rjval) + if (index(halofld,'L1') > 0) then + ! logical math doesn't work this way, force to correct answer + cichk = aichk ! p5 * (mod(nint(cidata_bas(i,j,k1,k2,iblock)),2) + rsign * mod(nint(rival),2)) + cjchk = ajchk ! p5 * (mod(nint(cidata_bas(i,j,k1,k2,iblock)),2) + rsign * mod(nint(rjval),2)) + else + cichk = p5 * (cidata_bas(i,j,k1,k2,iblock) + rsign * rival) + cjchk = p5 * (cjdata_bas(i,j,k1,k2,iblock) + rsign * rjval) + endif + else ! standard tripole fold cichk = rsign * rival cjchk = rsign * rjval endif -! if (testcnt == 6 .and. j == 61 .and. i < 3) then -! if (testcnt == 186 .and. j == 61 .and. i<4) then -! if (testcnt == 13 .and. j > 61 .and. (i < 3 .or. i > 89)) then -! if (testcnt == 5 .and. j >= 61 .and. (i < 3 .or. i > 90)) then -! write(100+my_task,'(a,5i6,2l3,f6.2,i6)') 'tcx1 ',i,j,iblock,itrip,jtrip, & -! tripole_average,tripole_pole,rsign,this_block%i_glob(i) -! write(100+my_task,'(a,4f12.2)') 'tcx2 ',cidata_std(i,j,k1,k2,iblock),rival,cichk,aichk -! write(100+my_task,'(a,4f12.2)') 'tcx3 ',cjdata_std(i,j,k1,k2,iblock),rjval,cjchk,ajchk -! endif endif ! tripole or tripoleT + endif if (index(halofld,'I4') > 0) then @@ -695,16 +713,16 @@ program halochk cjchk = real(nint(cjchk),kind=dbl_kind) endif - if (index(halofld,'L1') > 0 .and. .not.spvalL1) then + if (index(halofld,'L1') > 0) then if (cichk == dhalofillval .or. cichk == fillval) then cichk = c0 else - cichk = c1 + cichk = mod(nint(cichk),2) endif if (cjchk == dhalofillval .or. cjchk == fillval) then cjchk = c1 else - cjchk = c0 + cjchk = mod(nint(cjchk),2) endif endif @@ -719,6 +737,7 @@ program halochk enddo ! j enddo ! iblock + enddo ! maxfills enddo ! maxtypes enddo ! maxlocs enddo ! maxtests @@ -746,10 +765,10 @@ program halochk do n = 1,tottest if (errorflag(n) == passflag) then tpcnt = tpcnt + 1 - write(6,*) 'PASS ',trim(teststring(n)),ptcnt(n),failcnt(n) + write(6,'(2a,2i8)') 'PASS ',trim(teststring(n)),ptcnt(n),failcnt(n) else tfcnt = tfcnt + 1 - write(6,*) 'FAIL ',trim(teststring(n)),ptcnt(n),failcnt(n) + write(6,'(2a,2i8)') 'FAIL ',trim(teststring(n)),ptcnt(n),failcnt(n) endif enddo write(6,*) ' ' @@ -793,8 +812,10 @@ subroutine chkresults(a1,r1,errorflag,testcnt,failcnt,i,j,k1,k2,iblock,first_cal character(len=*) , parameter :: subname='(chkresults)' if (a1 /= r1 .or. print_always) then - errorflag = failflag - failcnt = failcnt + 1 + if (a1 /= r1) then + errorflag = failflag + failcnt = failcnt + 1 + endif if (first_call) then write(100+my_task,*) ' ' write(100+my_task,'(a,i4,2a)') '------- TEST = ',testcnt,' ',trim(teststring) diff --git a/cicecore/shared/ice_arrays_column.F90 b/cicecore/shared/ice_arrays_column.F90 index 38f3ee0f7..8dee4aef3 100644 --- a/cicecore/shared/ice_arrays_column.F90 +++ b/cicecore/shared/ice_arrays_column.F90 @@ -9,6 +9,7 @@ module ice_arrays_column use ice_kinds_mod + use ice_constants, only : c0 use ice_fileunits, only: nu_diag use ice_blocks, only: nx_block, ny_block use ice_domain_size, only: max_blocks, ncat, nilyr, nslyr, & @@ -25,8 +26,7 @@ module ice_arrays_column ! icepack_atmo.F90 ! Cdn variables on the T-grid - real (kind=dbl_kind), public, & - dimension (:,:,:), allocatable :: & + real (kind=dbl_kind), public, dimension (:,:,:), allocatable :: & Cdn_atm , & ! atm drag coefficient Cdn_ocn , & ! ocn drag coefficient ! form drag @@ -64,16 +64,17 @@ module ice_arrays_column ! icepack_itd.F90 real (kind=dbl_kind), public, allocatable :: & - hin_max(:) ! category limits (m) + hin_max(:) ! category limits (m) - character (len=35), public, allocatable :: c_hi_range(:) + character (len=35), public, allocatable :: & + c_hi_range(:)! string for history output ! icepack_snow.F90 real (kind=dbl_kind), public, dimension (:,:,:), allocatable :: & meltsliq ! snow melt mass (kg/m^2/step-->kg/m^2/day) real (kind=dbl_kind), public, dimension (:,:,:,:), allocatable :: & - meltsliqn ! snow melt mass in category n (kg/m^2) + meltsliqn ! snow melt mass in category n (kg/m^2) ! icepack_meltpond_lvl.F90 real (kind=dbl_kind), public, dimension (:,:,:,:), allocatable :: & @@ -83,10 +84,10 @@ module ice_arrays_column ! icepack_shortwave.F90 ! category albedos real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: & - alvdrn , & ! visible direct albedo (fraction) - alidrn , & ! near-ir direct albedo (fraction) - alvdfn , & ! visible diffuse albedo (fraction) - alidfn ! near-ir diffuse albedo (fraction) + alvdrn, & ! visible direct albedo (fraction) + alidrn, & ! near-ir direct albedo (fraction) + alvdfn, & ! visible diffuse albedo (fraction) + alidfn ! near-ir diffuse albedo (fraction) ! albedo components for history real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: & @@ -100,14 +101,14 @@ module ice_arrays_column ! shortwave components real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable, public :: & - Iswabsn ! SW radiation absorbed in ice layers (W m-2) + Iswabsn ! SW radiation absorbed in ice layers (W m-2) real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable, public :: & - Sswabsn ! SW radiation absorbed in snow layers (W m-2) + Sswabsn ! SW radiation absorbed in snow layers (W m-2) real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: & - fswsfcn , & ! SW absorbed at ice/snow surface (W m-2) - fswthrun , & ! SW through ice to ocean (W/m^2) + fswsfcn , & ! SW absorbed at ice/snow surface (W m-2) + fswthrun , & ! SW through ice to ocean (W/m^2) fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) @@ -119,7 +120,7 @@ module ice_arrays_column fswintn ! SW absorbed in ice interior, below surface (W m-2) real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable, public :: & - fswpenln ! visible SW entering ice layers (W m-2) + fswpenln ! visible SW entering ice layers (W m-2) ! biogeochemistry components @@ -348,6 +349,71 @@ subroutine alloc_arrays_column stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of Memory1') + Cdn_atm = c0 + Cdn_ocn = c0 + hfreebd = c0 + hdraft = c0 + hridge = c0 + distrdg = c0 + hkeel = c0 + dkeel = c0 + lfloe = c0 + dfloe = c0 + Cdn_atm_skin = c0 + Cdn_atm_floe = c0 + Cdn_atm_pond = c0 + Cdn_atm_rdg = c0 + Cdn_ocn_skin = c0 + Cdn_ocn_floe = c0 + Cdn_ocn_keel = c0 + Cdn_atm_ratio = c0 + grow_net = c0 + PP_net = c0 + hbri = c0 + chl_net = c0 + NO_net = c0 + upNO = c0 + upNH = c0 + meltsliq = c0 + meltsliqn = c0 + dhsn = c0 + ffracn = c0 + alvdrn = c0 + alidrn = c0 + alvdfn = c0 + alidfn = c0 + albicen = c0 + albsnon = c0 + albpndn = c0 + apeffn = c0 + snowfracn = c0 + fswsfcn = c0 + fswthrun = c0 + fswthrun_vdr = c0 + fswthrun_vdf = c0 + fswthrun_idr = c0 + fswthrun_idf = c0 + fswthrun_uvrdr= c0 + fswthrun_uvrdf= c0 + fswthrun_pardr= c0 + fswthrun_pardf= c0 + fswintn = c0 + first_ice_real= c0 + first_ice = .false. + dhbr_top = c0 + dhbr_bot = c0 + darcy_V = c0 + sice_rho = c0 + Iswabsn = c0 + Sswabsn = c0 + fswpenln = c0 + Zoo = c0 + zfswin = c0 + iDi = c0 + iki = c0 + bphi = c0 + bTiz = c0 + allocate( & ocean_bio (nx_block,ny_block,max_nbtrcr,max_blocks), & ! contains all the ocean bgc tracer concentrations fbio_snoice (nx_block,ny_block,max_nbtrcr,max_blocks), & ! fluxes from snow to ice @@ -359,6 +425,14 @@ subroutine alloc_arrays_column stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of Memory2') + ocean_bio = c0 + fbio_snoice = c0 + fbio_atmice = c0 + ocean_bio_all= c0 + ice_bio_net = c0 + snow_bio_net = c0 + algal_peak = 0 + allocate( & hin_max(0:ncat) , & ! category limits (m) c_hi_range(ncat) , & ! @@ -370,6 +444,14 @@ subroutine alloc_arrays_column stat=ierr) if (ierr/=0) call abort_ice(subname//' Out of Memory3') + hin_max = c0 + c_hi_range = '' + bgrid = c0 + igrid = c0 + cgrid = c0 + icgrid = c0 + swgrid = c0 + ! floe size distribution allocate( & floe_rad_l (nfsd) , & ! fsd size lower bound in m (radius) @@ -388,6 +470,20 @@ subroutine alloc_arrays_column stat=ierr) if (ierr/=0) call abort_ice(subname//' Out of Memory5') + floe_rad_l = c0 + floe_rad_c = c0 + floe_binwidth = c0 + c_fsd_range = '' + wavefreq = c0 + dwavefreq = c0 + wave_sig_ht = c0 + wave_spectrum = c0 + d_afsd_newi = c0 + d_afsd_latg = c0 + d_afsd_latm = c0 + d_afsd_wave = c0 + d_afsd_weld = c0 + end subroutine alloc_arrays_column !=======================================================================