diff --git a/components/elm/src/cpl/lnd_comp_mct.F90 b/components/elm/src/cpl/lnd_comp_mct.F90 index 47853d4c1cd8..9035afe2a079 100644 --- a/components/elm/src/cpl/lnd_comp_mct.F90 +++ b/components/elm/src/cpl/lnd_comp_mct.F90 @@ -406,9 +406,11 @@ subroutine lnd_init_mct( EClock, cdata_l, x2l_l, l2x_l, NLFilename ) ! write out the mesh file to disk, in parallel outfile = 'wholeLnd.h5m'//C_NULL_CHAR wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR - ierr = iMOAB_WriteMesh(mlnid, outfile, wopts) - if (ierr > 0 ) & - call endrun('Error: fail to write the land mesh file') + if (mlnid >= 0) then + ierr = iMOAB_WriteMesh(mlnid, outfile, wopts) + if (ierr > 0 ) & + call endrun('Error: fail to write the land mesh file') + endif #endif end subroutine lnd_init_mct @@ -1076,6 +1078,7 @@ subroutine lnd_export_moab(EClock, bounds, lnd2atm_vars, lnd2glc_vars) l2x_lm(i,index_l2x_Sl_tref) = lnd2atm_vars%t_ref2m_grc(g) l2x_lm(i,index_l2x_Sl_qref) = lnd2atm_vars%q_ref2m_grc(g) l2x_lm(i,index_l2x_Sl_u10) = lnd2atm_vars%u_ref10m_grc(g) + l2x_lm(i,index_l2x_Sl_u10withgusts)=lnd2atm_vars%u_ref10m_with_gusts_grc(g) ! see commit 5813d4103 l2x_lm(i,index_l2x_Fall_taux) = -lnd2atm_vars%taux_grc(g) l2x_lm(i,index_l2x_Fall_tauy) = -lnd2atm_vars%tauy_grc(g) l2x_lm(i,index_l2x_Fall_lat) = -lnd2atm_vars%eflx_lh_tot_grc(g) @@ -1339,7 +1342,7 @@ subroutine lnd_import_moab(EClock, bounds, atm2lnd_vars, glc2lnd_vars) #endif tagname=trim(seq_flds_x2l_fields)//C_NULL_CHAR ent_type = 0 ! vertices - ierr = iMOAB_GetDoubleTagStorage ( mlnid, tagname, totalmblsimp , ent_type, x2l_lm(1,1) ) + ierr = iMOAB_GetDoubleTagStorage ( mlnid, tagname, totalmblsimp , ent_type, x2l_lm ) if ( ierr > 0) then call endrun('Error: fail to get seq_flds_x2l_fields for land moab instance on component') endif diff --git a/components/mosart/src/cpl/rof_comp_mct.F90 b/components/mosart/src/cpl/rof_comp_mct.F90 index 9f0aa9fa405b..84eebe40c483 100644 --- a/components/mosart/src/cpl/rof_comp_mct.F90 +++ b/components/mosart/src/cpl/rof_comp_mct.F90 @@ -388,9 +388,11 @@ subroutine rof_init_mct( EClock, cdata_r, x2r_r, r2x_r, NLFilename) ! write out the mesh file to disk, in parallel outfile = 'wholeRof.h5m'//C_NULL_CHAR wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR - ierr = iMOAB_WriteMesh(mrofid, outfile, wopts) - if (ierr > 0 ) & - call shr_sys_abort( sub//' Error: fail to write the moab runoff mesh file') + if (mrofid >= 0) then + ierr = iMOAB_WriteMesh(mrofid, outfile, wopts) + if (ierr > 0 ) & + call shr_sys_abort( sub//' Error: fail to write the moab runoff mesh file') + endif #endif end subroutine rof_init_mct @@ -1260,7 +1262,7 @@ subroutine rof_import_moab( EClock ) ! ! LOCAL VARIABLES - integer :: n2, n, nt, begr, endr, nliq, nfrz + integer :: n2, n, nt, begr, endr, nliq, nfrz, nmud, nsan real(R8) :: tmp1, tmp2 real(R8) :: shum character(CXX) :: tagname ! @@ -1287,7 +1289,7 @@ subroutine rof_import_moab( EClock ) ! populate the array x2r_rm with data from MOAB tags tagname=trim(seq_flds_x2r_fields)//C_NULL_CHAR ent_type = 0 ! vertices, point cloud - ierr = iMOAB_GetDoubleTagStorage ( mrofid, tagname, totalmbls_r , ent_type, x2r_rm(1,1) ) + ierr = iMOAB_GetDoubleTagStorage ( mrofid, tagname, totalmbls_r , ent_type, x2r_rm ) if ( ierr > 0) then call shr_sys_abort(sub//'Error: fail to get seq_flds_a2x_fields for atm physgrid moab mesh') endif @@ -1295,6 +1297,8 @@ subroutine rof_import_moab( EClock ) ! Note that ***runin are fluxes nliq = 0 nfrz = 0 + nmud = 0 + nsan = 0 do nt = 1,nt_rtm if (trim(rtm_tracers(nt)) == 'LIQ') then nliq = nt @@ -1302,9 +1306,27 @@ subroutine rof_import_moab( EClock ) if (trim(rtm_tracers(nt)) == 'ICE') then nfrz = nt endif + if (trim(rtm_tracers(nt)) == 'MUD') then + nmud = nt + endif + if (trim(rtm_tracers(nt)) == 'SAN') then + nsan = nt + endif enddo - if (nliq == 0 .or. nfrz == 0) then - write(iulog,*) trim(sub),': ERROR in rtm_tracers LIQ ICE ',nliq,nfrz,rtm_tracers + if (nliq == 0) then + write(iulog,*) trim(sub),': ERROR in rtm_tracers LIQ',nliq,rtm_tracers + call shr_sys_abort() + endif + if (nfrz == 0) then + write(iulog,*) trim(sub),': ERROR in rtm_tracers ICE',nfrz,rtm_tracers + call shr_sys_abort() + endif + if (nmud == 0) then + write(iulog,*) trim(sub),': ERROR in rtm_tracers MUD',nmud,rtm_tracers + call shr_sys_abort() + endif + if (nsan == 0) then + write(iulog,*) trim(sub),': ERROR in rtm_tracers SAN',nsan,rtm_tracers call shr_sys_abort() endif @@ -1351,8 +1373,25 @@ subroutine rof_import_moab( EClock ) THeat%forc_vp(n) = shum * THeat%forc_pbot(n) / (0.622_r8 + 0.378_r8 * shum) THeat%coszen(n) = x2r_rm(n2,index_x2r_coszen_str) end if + + + rtmCTL%qsur(n,nmud) = 0.0_r8 + rtmCTL%qsur(n,nsan) = 0.0_r8 + + if (index_x2r_Flrl_inundinf > 0) then + rtmCTL%inundinf(n) = x2r_rm(n2,index_x2r_Flrl_inundinf) * (rtmCTL%area(n)*0.001_r8) + endif + enddo + if(sediflag) then + do n = begr,endr + n2 = n - begr + 1 + rtmCTL%qsur(n,nmud) = x2r_rm(n2,index_x2r_Flrl_rofmud) * (rtmCTL%area(n)) ! kg/m2/s --> kg/s for sediment + rtmCTL%qsur(n,nsan) = 0.0_r8 + enddo + end if + end subroutine rof_import_moab diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F index 8b8dac047234..c9465ca51983 100644 --- a/components/mpas-ocean/driver/ocn_comp_mct.F +++ b/components/mpas-ocean/driver/ocn_comp_mct.F @@ -3692,7 +3692,7 @@ subroutine ocn_import_moab(Eclock, errorCode)!{{{ ent_type = 1 ! cells ! get all tags in one method tagname = trim(seq_flds_x2o_fields)//C_NULL_CHAR - ierr = iMOAB_GetDoubleTagStorage ( MPOID, tagname, totalmbls_r , ent_type, x2o_om(1, 1) ) + ierr = iMOAB_GetDoubleTagStorage ( MPOID, tagname, totalmbls_r , ent_type, x2o_om ) if ( ierr /= 0 ) then write(ocnLogUnit,*) 'Fail to get MOAB fields ' endif diff --git a/components/mpas-seaice/driver/ice_comp_mct.F b/components/mpas-seaice/driver/ice_comp_mct.F index c42efa0a634e..2576624cac61 100644 --- a/components/mpas-seaice/driver/ice_comp_mct.F +++ b/components/mpas-seaice/driver/ice_comp_mct.F @@ -3550,6 +3550,7 @@ subroutine ice_export_moab(EClock) i2x_im(n, index_i2x_Si_t) = Tsrf i2x_im(n, index_i2x_Si_bpress) = basalPressure i2x_im(n, index_i2x_Si_u10) = atmosReferenceSpeed10m(i) + i2x_im(n, index_i2x_Si_u10withgusts) = atmosReferenceSpeed10m(i) ! see commit 5813d4103 i2x_im(n, index_i2x_Si_tref) = atmosReferenceTemperature2m(i) i2x_im(n, index_i2x_Si_qref) = atmosReferenceHumidity2m(i) i2x_im(n, index_i2x_Si_snowh) = snowVolumeCell(i) / ailohi @@ -3887,7 +3888,7 @@ subroutine ice_import_moab(Eclock)!{{{ ent_type = 1 ! cells ! set all tags in one method tagname = trim(seq_flds_x2i_fields)//C_NULL_CHAR - ierr = iMOAB_GetDoubleTagStorage ( MPSIID, tagname, totalmblr , ent_type, x2i_im(1, 1) ) + ierr = iMOAB_GetDoubleTagStorage ( MPSIID, tagname, totalmblr , ent_type, x2i_im ) if ( ierr /= 0 ) then write(iceLogUnit,*) 'Fail to get seq_flds_x2i_fields ' endif diff --git a/driver-moab/cime_config/buildnml b/driver-moab/cime_config/buildnml index 2abe08114108..ce3ed839635d 100755 --- a/driver-moab/cime_config/buildnml +++ b/driver-moab/cime_config/buildnml @@ -41,6 +41,7 @@ def _create_drv_namelists(case, infile, confdir, nmlgen, files): config['CPL_EPBAL'] = case.get_value('CPL_EPBAL') config['FLDS_WISO'] = case.get_value('FLDS_WISO') config['FLDS_POLAR'] = case.get_value('FLDS_POLAR') + config['FLDS_TF'] = case.get_value('FLDS_TF') config['BUDGETS'] = case.get_value('BUDGETS') config['MACH'] = case.get_value('MACH') config['MPILIB'] = case.get_value('MPILIB') @@ -69,6 +70,20 @@ def _create_drv_namelists(case, infile, confdir, nmlgen, files): elif case.get_value('RUN_TYPE') == 'branch': config['run_type'] = 'branch' + # --------------------------------------------------- + # set wave coupling settings based on compset: + # --------------------------------------------------- + if case.get_value('COMP_WAV') == 'ww3': + config['WAVSPEC'] = case.get_value('WAV_SPEC') + if case.get_value('COMP_OCN') == 'mpaso': + config['WAV_OCN_COUP'] = 'twoway' + elif case.get_value('COMP_OCN') == 'docn': + config['WAV_OCN_COUP'] = 'oneway' + elif case.get_value('COMP_WAV') == 'dwav': + config['WAVSPEC'] = 'sp36x36' + else: + config['WAVSPEC'] = 'none' + #---------------------------------------------------- # Initialize namelist defaults #---------------------------------------------------- diff --git a/driver-moab/cime_config/namelist_definition_drv.xml b/driver-moab/cime_config/namelist_definition_drv.xml index 8b633b1100be..12117982b466 100644 --- a/driver-moab/cime_config/namelist_definition_drv.xml +++ b/driver-moab/cime_config/namelist_definition_drv.xml @@ -149,6 +149,18 @@ + + logical + seq_flds + seq_cplflds_inparm + + If set to .true. thermal forcing fields will be passed from the ocean to the coupler. + + + $FLDS_TF + + + logical seq_flds @@ -302,6 +314,18 @@ .false. + + + char + seq_flds + seq_cplflds_inparm + One- or Two-way coupling between Wave and Ocn. + + none + oneway + twoway + + diff --git a/driver-moab/main/cime_comp_mod.F90 b/driver-moab/main/cime_comp_mod.F90 index ad0edae64edb..e95cddda8c22 100644 --- a/driver-moab/main/cime_comp_mod.F90 +++ b/driver-moab/main/cime_comp_mod.F90 @@ -452,6 +452,8 @@ module cime_comp_mod logical :: lnd_c2_glc ! .true. => lnd to glc coupling on logical :: ocn_c2_atm ! .true. => ocn to atm coupling on logical :: ocn_c2_ice ! .true. => ocn to ice coupling on + logical :: ocn_c2_glctf ! .true. => ocn to glc thermal forcing coupling on + integer :: glc_nzoc ! number of z-levels for ocn/glc TF coupling logical :: ocn_c2_glcshelf ! .true. => ocn to glc ice shelf coupling on logical :: ocn_c2_wav ! .true. => ocn to wav coupling on logical :: ocn_c2_rof ! .true. => ocn to rof coupling on @@ -569,7 +571,7 @@ module cime_comp_mod 'Sa_u:Sa_v' character(CL) :: hist_a2x3hr_flds = & - 'Sa_z:Sa_u:Sa_v:Sa_tbot:Sa_ptem:Sa_shum:Sa_dens:Sa_pbot:Sa_pslv:Faxa_lwdn:& + 'Sa_z:Sa_topo:Sa_u:Sa_v:Sa_tbot:Sa_ptem:Sa_shum:Sa_dens:Sa_pbot:Sa_pslv:Faxa_lwdn:& &Faxa_rainc:Faxa_rainl:Faxa_snowc:Faxa_snowl:& &Faxa_swndr:Faxa_swvdr:Faxa_swndf:Faxa_swvdf:& &Sa_co2diag:Sa_co2prog' @@ -1705,6 +1707,8 @@ subroutine cime_init() ocn_prognostic=ocn_prognostic, & ocnrof_prognostic=ocnrof_prognostic, & ocn_c2_glcshelf=ocn_c2_glcshelf, & + ocn_c2_glctf=ocn_c2_glctf, & + glc_nzoc=glc_nzoc, & glc_prognostic=glc_prognostic, & rof_prognostic=rof_prognostic, & rofocn_prognostic=rofocn_prognostic, & @@ -1801,6 +1805,7 @@ subroutine cime_init() if (atm_prognostic) ocn_c2_atm = .true. if (atm_present ) ocn_c2_atm = .true. ! needed for aoflux calc if aoflux=atm if (ice_prognostic) ocn_c2_ice = .true. + if (glc_prognostic .and. (glc_nzoc > 0)) ocn_c2_glctf = .true. if (wav_prognostic) ocn_c2_wav = .true. if (rofocn_prognostic) ocn_c2_rof = .true. @@ -2539,7 +2544,7 @@ subroutine cime_init() call shr_sys_flush(logunit) end if call t_startf('CPL:seq_rest_read-moab') - call seq_rest_mb_read(rest_file, infodata, samegrid_al) + call seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) call t_stopf('CPL:seq_rest_read-moab') #ifdef MOABDEBUG call write_moab_state(.false.) @@ -3529,7 +3534,7 @@ subroutine cime_run() call shr_sys_flush(logunit) end if call t_startf('CPL:seq_rest_read-moab') - call seq_rest_mb_read(drv_resume_file, infodata, samegrid_al) + call seq_rest_mb_read(drv_resume_file, infodata, samegrid_al, samegrid_lr) call t_stopf('CPL:seq_rest_read-moab') end if ! Clear the resume file so we don't try to read it again @@ -5388,7 +5393,7 @@ subroutine cime_run_write_restart(drv_pause, write_restart, drv_resume_file) call t_startf('CPL:seq_rest_mb_write') call seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & atm, lnd, ice, ocn, rof, glc, wav, esp, iac, & - trim(cpl_inst_tag), samegrid_al, drv_moab_resume_file) + trim(cpl_inst_tag), samegrid_al, samegrid_lr, drv_moab_resume_file) call t_stopf('CPL:seq_rest_mb_write') if (iamroot_CPLID) then diff --git a/driver-moab/main/component_mod.F90 b/driver-moab/main/component_mod.F90 index c250634cad6c..8d1bd2081e23 100644 --- a/driver-moab/main/component_mod.F90 +++ b/driver-moab/main/component_mod.F90 @@ -470,6 +470,7 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, use iso_c_binding ! character(1024) :: domain_file ! file containing domain info (set my input) use seq_comm_mct, only: mboxid ! iMOAB id for MPAS ocean migrated mesh to coupler pes + use seq_comm_mct, only: mblxid ! iMOAB id for lnd migrated mesh to coupler pes use seq_comm_mct, only: mbaxid ! iMOAB id for atm migrated mesh to coupler pes use seq_comm_mct, only: mbrxid ! iMOAB id for rof migrated mesh to coupler pes #endif @@ -551,18 +552,6 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, ! project now aream on ocean (from atm) #endif call seq_map_map(mapper_Fa2o, av_s=dom_s%data, av_d=dom_d%data, fldlist='aream') - -#ifdef HAVE_MOAB -#ifdef MOABDEBUG - ierr = iMOAB_WriteMesh(mboxid, trim('recMeshOcnWithArea.h5m'//C_NULL_CHAR), & - trim(';PARALLEL=WRITE_PART'//C_NULL_CHAR)) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing ocean mesh coupler ' - call shr_sys_abort(subname//' ERROR in writing ocean mesh coupler ') - endif -#endif -#endif - else gsmap_s => component_get_gsmap_cx(ocn(1)) ! gsmap_ox @@ -644,6 +633,40 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, endif endif +#ifdef MOABDEBUG + if(mbaxid >=0 ) then + ierr = iMOAB_WriteMesh(mbaxid, trim('cplAtmWithAream.h5m'//C_NULL_CHAR), & + trim(';PARALLEL=WRITE_PART'//C_NULL_CHAR)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing atm mesh coupler ' + call shr_sys_abort(subname//' ERROR in writing atm mesh coupler ') + endif + endif + if(mblxid >=0 ) then + ierr = iMOAB_WriteMesh(mblxid, trim('cplLndWithAream.h5m'//C_NULL_CHAR), & + trim(';PARALLEL=WRITE_PART'//C_NULL_CHAR)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing lnd mesh coupler ' + call shr_sys_abort(subname//' ERROR in writing lnd mesh coupler ') + endif + endif + if(mboxid >=0 ) then + ierr = iMOAB_WriteMesh(mboxid, trim('cplOcnWithAream.h5m'//C_NULL_CHAR), & + trim(';PARALLEL=WRITE_PART'//C_NULL_CHAR)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing ocn mesh coupler ' + call shr_sys_abort(subname//' ERROR in writing ocn mesh coupler ') + endif + endif + if(mbrxid >=0 ) then + ierr = iMOAB_WriteMesh(mbrxid, trim('cplRofWithAream.h5m'//C_NULL_CHAR), & + trim(';PARALLEL=WRITE_PART'//C_NULL_CHAR)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing rof mesh coupler ' + call shr_sys_abort(subname//' ERROR in writing rof mesh coupler ') + endif + endif +#endif end subroutine component_init_aream @@ -734,18 +757,22 @@ subroutine component_init_areacor_moab (comp, mbccid, mbcxid, seq_flds_c2x_fluxe integer :: mpi_tag character(*), parameter :: subname = '(component_init_areacor_moab)' character(CXX) :: tagname + character(CXX) :: comment integer :: tagtype, numco, tagindex, lsize, i, j, arrsize, ierr, nfields real (kind=r8) , allocatable :: areas (:,:), factors(:,:), vals(:,:) ! 2 tags values, area, aream, - real (kind=r8) :: rarea, raream, rmask, fact + real (kind=r8) :: rarea, raream, rmask, fact, rmin1, rmax1, rmin, rmax integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) type(mct_list) :: temp_list ! used to count number of fields + integer :: mpicom + logical :: iamroot + character(len=*),parameter :: F0R = "(2A,2g23.15,A )" !--------------------------------------------------------------- if (comp(1)%iamin_cplcompid) then tagname='aream'//C_NULL_CHAR ! bring on the comp side the aream from maps ! (it is either computed by mapping routine or read from mapping files) - call component_exch_moab(comp(1), mbcxid, mbccid, 1, tagname) + call component_exch_moab(comp(1), mbcxid, mbccid, 1, tagname, context_exch='aream') ! For only component pes if (comp(1)%iamin_compid) then @@ -795,6 +822,23 @@ subroutine component_init_areacor_moab (comp, mbccid, mbcxid, seq_flds_c2x_fluxe call shr_sys_abort(subname//' cannot set correction area factors ') endif +! print the computed values for area factors + rmin1 = minval(factors(:,1)) + rmax1 = maxval(factors(:,1)) + mpicom = comp(1)%mpicom_compid + iamroot= comp(1)%iamroot_compid + call shr_mpi_min(rmin1,rmin,mpicom) + call shr_mpi_max(rmax1,rmax,mpicom) + comment = 'areafact_'//trim(comp(1)%name) + if (iamroot) write(logunit,F0R) trim(subname),' : min/max mdl2drv ',rmin,rmax,trim(comment) + + rmin1 = minval(factors(:,2)) + rmax1 = maxval(factors(:,2)) + call shr_mpi_min(rmin1,rmin,mpicom) + call shr_mpi_max(rmax1,rmax,mpicom) + if (iamroot) write(logunit,F0R) trim(subname),' : min/max drv2mdl ',rmin,rmax,trim(comment) + if (iamroot) call shr_sys_flush(logunit) + ! Area correct component initialization output fields ! need to multiply fluxes (correct them) with mdl2drv (factors(i,1)) ! so get all fluxes (tags) multiply with factor(i,1), according to mask @@ -836,7 +880,7 @@ subroutine component_init_areacor_moab (comp, mbccid, mbcxid, seq_flds_c2x_fluxe endif ! send data to coupler exchange ? everything, not only fluxes ? - call component_exch_moab(comp(1), mbccid, mbcxid, 0, seq_flds_c2x_fields) + call component_exch_moab(comp(1), mbccid, mbcxid, 0, seq_flds_c2x_fields, context_exch='areacor') endif end subroutine component_init_areacor_moab diff --git a/driver-moab/main/component_type_mod.F90 b/driver-moab/main/component_type_mod.F90 index ebb3ad8f58e3..e42acf207d38 100644 --- a/driver-moab/main/component_type_mod.F90 +++ b/driver-moab/main/component_type_mod.F90 @@ -421,6 +421,7 @@ subroutine compare_mct_av_moab_tag(comp, attrVect, mct_field, appId, tagname, en use shr_mpi_mod, only: shr_mpi_sum use shr_kind_mod, only: CXX => shr_kind_CXX + use shr_kind_mod , only: IN=>SHR_KIND_IN use seq_comm_mct , only : CPLID, seq_comm_iamroot use seq_comm_mct, only: seq_comm_setptrs use iMOAB, only : iMOAB_DefineTagStorage, iMOAB_GetDoubleTagStorage, & @@ -431,20 +432,22 @@ subroutine compare_mct_av_moab_tag(comp, attrVect, mct_field, appId, tagname, en type(component_type), intent(in) :: comp integer , intent(in) :: appId, ent_type type(mct_aVect) , intent(in), pointer :: attrVect + type(mct_gsmap), pointer :: gsmap character(*) , intent(in) :: mct_field character(*) , intent(in) :: tagname + integer(IN), pointer :: dof(:) real(r8) , intent(out) :: difference logical , intent(in) :: first_time real(r8) :: differenceg ! global, reduced diff - type(mct_ggrid), pointer :: dom - integer :: kgg, mbSize, nloc, index_avfield + !type(mct_ggrid), pointer :: dom + integer :: kgg, mbSize, nloc, index_avfield, my_task ! moab integer :: tagtype, numco, tagindex, ierr character(CXX) :: tagname_mct - integer , allocatable :: GlobalIds(:) ! used for setting values associated with ids + !integer , allocatable :: GlobalIds(:) ! used for setting values associated with ids real(r8) , allocatable :: values(:), mct_values(:) integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) @@ -457,11 +460,15 @@ subroutine compare_mct_av_moab_tag(comp, attrVect, mct_field, appId, tagname, en call seq_comm_setptrs(CPLID, mpicom=mpicom) nloc = mct_avect_lsize(attrVect) - allocate(GlobalIds(nloc)) + !allocate(GlobalIds(nloc)) allocate(values(nloc)) - dom => component_get_dom_cx(comp) - kgg = mct_aVect_indexIA(dom%data ,"GlobGridNum" ,perrWith=subName) - GlobalIds = dom%data%iAttr(kgg,:) + !dom => component_get_dom_cx(comp) + !kgg = mct_aVect_indexIA(dom%data ,"GlobGridNum" ,perrWith=subName) + !GlobalIds = dom%data%iAttr(kgg,:) + gsmap => component_get_gsmap_cx(comp) ! gsmap_x + call mpi_comm_rank(mpicom,my_task,ierr) + ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid + call mct_gsMap_orderedPoints(gsmap, my_task, dof) index_avfield = mct_aVect_indexRA(attrVect,trim(mct_field)) values(:) = attrVect%rAttr(index_avfield,:) @@ -476,7 +483,7 @@ subroutine compare_mct_av_moab_tag(comp, attrVect, mct_field, appId, tagname, en if (ierr > 0 ) & call shr_sys_abort(subname//'Error: fail to define new tag for mct') endif - ierr = iMOAB_SetDoubleTagStorageWithGid ( appId, tagname_mct, nloc , ent_type, values, GlobalIds ) + ierr = iMOAB_SetDoubleTagStorageWithGid ( appId, tagname_mct, nloc , ent_type, values, dof ) if (ierr > 0 ) & call shr_sys_abort(subname//'Error: fail to set new tags') @@ -511,7 +518,8 @@ subroutine compare_mct_av_moab_tag(comp, attrVect, mct_field, appId, tagname, en print * , subname, trim(comp%ntype), ' on cpl, difference on tag ', trim(tagname), ' = ', difference !call shr_sys_abort(subname//'differences between mct and moab values') endif - deallocate(GlobalIds) + !deallocate(GlobalIds) + deallocate (dof) deallocate(values) deallocate(mct_values) diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 8975a9e9e62f..04557558ea1c 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -60,6 +60,7 @@ module cplcomp_exchange_mod private :: seq_mctext_gsmapCreate private :: seq_mctext_avCreate + private :: copy_aream_from_area !-------------------------------------------------------------------------- ! Public data !-------------------------------------------------------------------------- @@ -981,7 +982,34 @@ logical function seq_mctext_gsmapIdentical(gsmap1,gsmap2) end function seq_mctext_gsmapIdentical +subroutine copy_aream_from_area(mbappid) + + ! maybe we will move this from here + use iMOAB, only: iMOAB_GetDoubleTagStorage, iMOAB_SetDoubleTagStorage, iMOAB_GetMeshInfo + + integer , intent(in) :: mbappid + character(CXX) :: tagname + integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) + real(r8), allocatable :: tagValues(:) ! used for setting aream tags for atm domain read case + integer :: arrSize ! for the size of tagValues + integer :: ierr, ent_type + + ! copy aream from area + if (mbappid >= 0) then ! coupler atm procs + ierr = iMOAB_GetMeshInfo ( mbappid, nvert, nvise, nbl, nsurf, nvisBC ) + arrSize = nvise(1) ! cells + allocate(tagValues(arrSize)) + tagname = 'area'//C_NULL_CHAR + ent_type = 1 ! cells + ierr = iMOAB_GetDoubleTagStorage( mbappid, tagname, arrsize , ent_type, tagValues ) + tagname = 'aream'//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage( mbappid, tagname, arrsize , ent_type, tagValues ) + deallocate(tagValues) + endif + + return !======================================================================= + end subroutine copy_aream_from_area subroutine cplcomp_moab_Init(infodata,comp) @@ -999,6 +1027,7 @@ subroutine cplcomp_moab_Init(infodata,comp) ! type(seq_infodata_type) , intent(in) :: infodata type(component_type), intent(inout) :: comp + ! ! Local Variables ! @@ -1015,7 +1044,7 @@ subroutine cplcomp_moab_Init(infodata,comp) integer :: mpigrp_cplid ! coupler pes integer :: mpigrp_old ! component group pes integer :: ierr, context_id - character*200 :: appname, outfile, wopts, ropts + character*200 :: appname, outfile, wopts, ropts, infile character(CL) :: rtm_mesh, rof_domain character(CL) :: lnd_domain character(CL) :: ocn_domain @@ -1028,6 +1057,8 @@ subroutine cplcomp_moab_Init(infodata,comp) ! and atm spectral on coupler character(CXX) :: tagname integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) + real(r8), allocatable :: tagValues(:) ! used for setting aream tags for atm domain read case + integer :: arrSize ! for the size of tagValues ! type(mct_list) :: temp_list ! integer :: nfields, arrsize ! real(R8), allocatable, target :: values(:) @@ -1112,15 +1143,17 @@ subroutine cplcomp_moab_Init(infodata,comp) endif else ! data atm ! we need to read the atm mesh on coupler, from domain file - ierr = iMOAB_LoadMesh(mbaxid, trim(atm_mesh)//C_NULL_CHAR, & - "PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;REPARTITION;NO_CULLING", 0) + infile = trim(atm_mesh)//C_NULL_CHAR + ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;REPARTITION;NO_CULLING'//C_NULL_CHAR + if (seq_comm_iamroot(CPLID)) then + write(logunit,'(A)') subname//' loading atm domain mesh from file '//trim(atm_mesh) & + , ' with options ' // trim(ropts) + endif + ierr = iMOAB_LoadMesh(mbaxid, infile, ropts, 0) if ( ierr /= 0 ) then write(logunit,*) 'Failed to load atm domain mesh on coupler' call shr_sys_abort(subname//' ERROR Failed to load atm domain mesh on coupler ') endif - if (seq_comm_iamroot(CPLID)) then - write(logunit,'(A)') subname//' load atm domain mesh from file '//trim(atm_mesh) - endif ! right now, turn atm_pg_active to true atm_pg_active = .true. ! FIXME TODO ! need to add global id tag to the app, it will be used in restart @@ -1202,8 +1235,9 @@ subroutine cplcomp_moab_Init(infodata,comp) ! also, frac, area, masks has to come from atm mphaid, not from domain file reader ! this is hard to digest :( tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR - call component_exch_moab(comp, mphaid, mbaxid, 0, tagname) - + call component_exch_moab(comp, mphaid, mbaxid, 0, tagname, context_exch='dom') + ! copy aream from area in case atm_mesh + call copy_aream_from_area(mbaxid) endif ! coupler pes #ifdef MOABDEBUG @@ -1284,11 +1318,13 @@ subroutine cplcomp_moab_Init(infodata,comp) endif else ! we need to read the ocean mesh on coupler, from domain file + infile = trim(ocn_domain)//C_NULL_CHAR + ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;NO_CULLING;REPARTITION'//C_NULL_CHAR if (seq_comm_iamroot(CPLID)) then - write(logunit,'(A)') subname//' loading ocn domain mesh from file '//trim(ocn_domain) + write(logunit,'(A)') subname//' loading ocn domain mesh from file '//trim(infile) & + , ' with options ' // trim(ropts) endif - ierr = iMOAB_LoadMesh(mboxid, trim(ocn_domain)//C_NULL_CHAR, & - "PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;NO_CULLING;REPARTITION", 0) + ierr = iMOAB_LoadMesh(mboxid, infile, ropts, 0) if ( ierr /= 0 ) then write(logunit,*) 'Failed to load ocean domain mesh on coupler' call shr_sys_abort(subname//' ERROR Failed to load ocean domain mesh on coupler ') @@ -1364,7 +1400,7 @@ subroutine cplcomp_moab_Init(infodata,comp) ! also, frac, area, masks has to come from ocean mpoid, not from domain file reader ! this is hard to digest :( tagname = 'area:frac:mask'//C_NULL_CHAR - call component_exch_moab(comp, mpoid, mboxid, 0, tagname) + call component_exch_moab(comp, mpoid, mboxid, 0, tagname, context_exch='afm') endif ! start copy @@ -1398,15 +1434,17 @@ subroutine cplcomp_moab_Init(infodata,comp) endif else ! we need to read the ocean mesh on coupler, from domain file - ierr = iMOAB_LoadMesh(mbofxid, trim(ocn_domain)//C_NULL_CHAR, & - "PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;NO_CULLING;REPARTITION", 0) + infile = trim(ocn_domain)//C_NULL_CHAR + ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;NO_CULLING;REPARTITION'//C_NULL_CHAR + if (seq_comm_iamroot(CPLID)) then + write(logunit,'(A)') subname//' load ocn domain mesh from file for second ocn instance '//trim(ocn_domain) & + , ' with options '//trim(ropts) + endif + ierr = iMOAB_LoadMesh(mbofxid, infile, ropts, 0) if ( ierr /= 0 ) then write(logunit,*) 'Failed to load second ocean domain mesh on coupler' call shr_sys_abort(subname//' ERROR Failed to load second ocean domain mesh on coupler ') endif - if (seq_comm_iamroot(CPLID)) then - write(logunit,'(A)') subname//' load ocn domain mesh from file for second ocn instance '//trim(ocn_domain) - endif ! need to add global id tag to the app, it will be used in restart tagtype = 0 ! dense, integer numco = 1 @@ -1441,6 +1479,7 @@ subroutine cplcomp_moab_Init(infodata,comp) endif ! end copy #ifdef MOABDEBUG + if (mbofxid >= 0) then outfile = 'recMeshOcnF.h5m'//C_NULL_CHAR wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! ! write out the mesh file to disk @@ -1449,6 +1488,7 @@ subroutine cplcomp_moab_Init(infodata,comp) write(logunit,*) subname,' error in writing ocean mesh coupler ' call shr_sys_abort(subname//' ERROR in writing ocean mesh coupler ') endif + endif #endif endif ! end ocean @@ -1475,7 +1515,7 @@ subroutine cplcomp_moab_Init(infodata,comp) outfile = trim(lnd_domain)//C_NULL_CHAR nghlay = 0 ! no ghost layers if (seq_comm_iamroot(CPLID) ) then - write(logunit, *) "load land domain file from file: ", trim(lnd_domain), & + write(logunit, *) "loading land domain file from file: ", trim(lnd_domain), & " with options: ", trim(ropts) endif ierr = iMOAB_LoadMesh(mblxid, outfile, ropts, nghlay) @@ -1483,9 +1523,6 @@ subroutine cplcomp_moab_Init(infodata,comp) write(logunit,*) subname,' error in reading land coupler mesh from ', trim(lnd_domain) call shr_sys_abort(subname//' ERROR in reading land coupler mesh') endif - if (seq_comm_iamroot(CPLID)) then - write(logunit,'(A)') subname//' load lnd domain mesh from file '//trim(lnd_domain) - endif ! need to add global id tag to the app, it will be used in restart tagtype = 0 ! dense, integer numco = 1 @@ -1540,7 +1577,7 @@ subroutine cplcomp_moab_Init(infodata,comp) endif tagname = 'lat:lon:area:frac:mask'//C_NULL_CHAR - call component_exch_moab(comp, mlnid, mblxid, 0, tagname) + call component_exch_moab(comp, mlnid, mblxid, 0, tagname, context_exch='dom') #ifdef MOABDEBUG outfile = 'recMeshLand.h5m'//C_NULL_CHAR @@ -1603,15 +1640,17 @@ subroutine cplcomp_moab_Init(infodata,comp) endif else ! we need to read the mesh ice (domain file) - ierr = iMOAB_LoadMesh(mbixid, trim(ice_domain)//C_NULL_CHAR, & - "PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;NO_CULLING;REPARTITION", 0) + ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;NO_CULLING;REPARTITION'//C_NULL_CHAR + infile = trim(ice_domain)//C_NULL_CHAR + if (seq_comm_iamroot(CPLID)) then + write(logunit,'(A)') subname//' loading ice domain mesh from file '//infile & + , ' with options '//trim(ropts) + endif + ierr = iMOAB_LoadMesh(mbixid, infile, ropts, 0) if ( ierr /= 0 ) then write(logunit,*) 'Failed to load ice domain mesh on coupler' call shr_sys_abort(subname//' ERROR Failed to load ice domain mesh on coupler ') endif - if (seq_comm_iamroot(CPLID)) then - write(logunit,'(A)') subname//' load ice domain mesh from file '//trim(ice_domain) - endif ! need to add global id tag to the app, it will be used in restart tagtype = 0 ! dense, integer numco = 1 @@ -1715,7 +1754,7 @@ subroutine cplcomp_moab_Init(infodata,comp) appname = "COUPLE_MROF"//C_NULL_CHAR ierr = iMOAB_RegisterApplication(trim(appname), mpicom_new, id_join, mbrxid) - ! load mesh from scrip file passed from river model + ! load mesh from scrip file passed from river model, if domain file is not available call seq_infodata_GetData(infodata,rof_mesh=rtm_mesh,rof_domain=rof_domain) if ( trim(rof_domain) == 'none' ) then outfile = trim(rtm_mesh)//C_NULL_CHAR @@ -1725,14 +1764,14 @@ subroutine cplcomp_moab_Init(infodata,comp) ropts = 'PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;REPARTITION'//C_NULL_CHAR endif nghlay = 0 ! no ghost layers - ierr = iMOAB_LoadMesh(mbrxid, outfile, ropts, nghlay) if (seq_comm_iamroot(CPLID)) then - write(logunit,'(A)') subname//' load rof from file '//trim(outfile) + write(logunit,'(A)') subname//' loading rof from file '//trim(outfile) & + , ' with options ', trim(ropts) endif + ierr = iMOAB_LoadMesh(mbrxid, outfile, ropts, nghlay) if ( ierr .ne. 0 ) then call shr_sys_abort( subname//' ERROR: cannot read rof mesh on coupler' ) end if - ! need to add global id tag to the app, it will be used in restart tagtype = 0 ! dense, integer numco = 1 @@ -1785,23 +1824,19 @@ subroutine cplcomp_moab_Init(infodata,comp) tagname = 'area:lon:lat:frac:mask'//C_NULL_CHAR call component_exch_moab(comp, mrofid, mbrxid, 0, tagname) - - ! if (mrofid .ge. 0) then ! we are on component rof pes - ! context_id = id_join - ! ierr = iMOAB_FreeSenderBuffers(mrofid, context_id) - ! if (ierr .ne. 0) then - ! write(logunit,*) subname,' error in freeing buffers ' - ! call shr_sys_abort(subname//' ERROR in freeing buffers ') - ! endif - ! endif + ! copy aream from area in all cases + ! initialize aream from area; it may have different values in the end, or reset again + call copy_aream_from_area(mbrxid) #ifdef MOABDEBUG - outfile = 'recMeshRof.h5m'//C_NULL_CHAR - wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! - ! write out the mesh file to disk - ierr = iMOAB_WriteMesh(mbrxid, trim(outfile), trim(wopts)) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing rof mesh on coupler ' - call shr_sys_abort(subname//' ERROR in writing rof mesh on coupler ') + if (mbrxid >= 0) then + outfile = 'recMeshRof.h5m'//C_NULL_CHAR + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! + ! write out the mesh file to disk + ierr = iMOAB_WriteMesh(mbrxid, trim(outfile), trim(wopts)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing rof mesh on coupler ' + call shr_sys_abort(subname//' ERROR in writing rof mesh on coupler ') + endif endif #endif endif ! end for rof coupler set up @@ -1811,7 +1846,7 @@ end subroutine cplcomp_moab_Init ! can exchange data between mesh in component and mesh on coupler. Either way. ! used in first hop of 2-hop - subroutine component_exch_moab(comp, mbAPPid1, mbAppid2, direction, fields ) + subroutine component_exch_moab(comp, mbAPPid1, mbAppid2, direction, fields, context_exch ) use iMOAB , only: iMOAB_SendElementTag, iMOAB_ReceiveElementTag, iMOAB_WriteMesh, iMOAB_FreeSenderBuffers use seq_comm_mct, only : num_moab_exports ! for debugging @@ -1825,6 +1860,7 @@ subroutine component_exch_moab(comp, mbAPPid1, mbAppid2, direction, fields ) ! direction 0 is from component to coupler; 1 is from coupler to component integer, intent(in) :: mbAPPid1, mbAppid2, direction character(CXX) , intent(in) :: fields + character(len=*) , intent(in), optional :: context_exch character(*), parameter :: subname = '(component_exch_moab)' integer :: id_join, source_id, target_id, ierr @@ -1854,8 +1890,25 @@ subroutine component_exch_moab(comp, mbAPPid1, mbAppid2, direction, fields ) if (comp%oneletterid == 'a' .and. direction .eq. 1 ) then target_id = target_id + 200 endif - if (mbAPPid1 .ge. 0) then ! send - + if (mbAPPid1 .ge. 0) then ! we are on the sending pes +#ifdef MOABDEBUG + if (direction .eq. 0 ) then + dir = 'c2x' + else + dir = 'x2c' + endif + write(lnum,"(I0.2)") num_moab_exports + if (present(context_exch)) then + outfile = comp%ntype//'_src_'//trim(context_exch)//'_'//trim(dir)//'_'//trim(lnum)//'.h5m'//C_NULL_CHAR + else + outfile = comp%ntype//'_src_'//trim(dir)//'_'//trim(lnum)//'.h5m'//C_NULL_CHAR + endif + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! + ierr = iMOAB_WriteMesh(mbAPPid1, trim(outfile), trim(wopts)) + if (ierr .ne. 0) then + call shr_sys_abort(subname//' cannot write file '// outfile) + endif +#endif ! basically, use the initial partitioning ierr = iMOAB_SendElementTag(mbAPPid1, tagName, mpicom_join, target_id) if (ierr .ne. 0) then @@ -1884,14 +1937,18 @@ subroutine component_exch_moab(comp, mbAPPid1, mbAppid2, direction, fields ) else dir = 'x2c' endif + write(lnum,"(I0.2)") num_moab_exports if (seq_comm_iamroot(CPLID) ) then - write(logunit,'(A)') subname//' '//comp%ntype//' called in direction '//trim(dir)//' for fields '//trim(tagname) + write(logunit,'(A)') subname//' '//comp%ntype//' at moab count '//trim(lnum)//' called in direction '//trim(dir)//' for fields '//trim(tagname) endif if (mbAPPid2 .ge. 0 ) then ! we are on receiving pes, for sure ! number_proj = number_proj+1 ! count the number of projections - write(lnum,"(I0.2)") num_moab_exports - outfile = comp%ntype//'_'//trim(dir)//'_'//trim(lnum)//'.h5m'//C_NULL_CHAR + if (present(context_exch)) then + outfile = comp%ntype//'_'//trim(context_exch)//'_'//trim(dir)//'_'//trim(lnum)//'.h5m'//C_NULL_CHAR + else + outfile = comp%ntype//'_'//trim(dir)//'_'//trim(lnum)//'.h5m'//C_NULL_CHAR + endif wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! ierr = iMOAB_WriteMesh(mbAPPid2, trim(outfile), trim(wopts)) if (ierr .ne. 0) then diff --git a/driver-moab/main/prep_atm_mod.F90 b/driver-moab/main/prep_atm_mod.F90 index 6e0d0ffef57d..7ea2560ab8df 100644 --- a/driver-moab/main/prep_atm_mod.F90 +++ b/driver-moab/main/prep_atm_mod.F90 @@ -127,7 +127,7 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at use iMOAB, only: iMOAB_ComputeMeshIntersectionOnSphere, iMOAB_RegisterApplication, & iMOAB_WriteMesh , iMOAB_ComputeCommGraph, iMOAB_ComputeScalarProjectionWeights, & - iMOAB_DefineTagStorage, iMOAB_ComputeCoverageMesh + iMOAB_DefineTagStorage, iMOAB_MigrateMapMesh !--------------------------------------------------------------- ! Description ! Initialize module attribute vectors and mappers @@ -275,6 +275,7 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at mapper_So2a%src_context = ocn(1)%cplcompid mapper_So2a%weight_identifier = wgtIdo2a mapper_So2a%mbname = 'mapper_So2a' + mapper_So2a%intx_context = idintx ! Since we are projecting fields from OCN to ATM-PHY grid, we need to define ! OCN o2x fields to ATM-PHY grid (or ATM-DYN (spectral) ) on coupler side @@ -294,14 +295,6 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at if (.not. samegrid_ao) then ! most cases - ! compute the OCN coverage mesh here on coupler pes - ! OCN mesh was redistributed to cover target (ATM) partition - ierr = iMOAB_ComputeCoverageMesh( mboxid, mbaxid, mbintxoa ) - if (ierr .ne. 0) then - write(logunit,*) subname,' cannot compute source OCN coverage mesh for ATM' - call shr_sys_abort(subname//' ERROR in computing OCN-ATM coverage') - endif - if (compute_maps_online_o2a) then if (iamroot_CPLID) then write(logunit,*) '....Computing weights' @@ -316,36 +309,29 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at endif #ifdef MOABDEBUG - wopts = C_NULL_CHAR - call shr_mpi_commrank( mpicom_CPLID, rank ) - if (rank .lt. 3) then - write(lnum,"(I0.2)")rank ! - outfile = 'intx_oa_'//trim(lnum)// '.h5m' // C_NULL_CHAR - ierr = iMOAB_WriteMesh(mbintxoa, outfile, wopts) ! write local intx file - if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing OCN-ATM intersection file ' - call shr_sys_abort(subname//' ERROR in writing OCN-ATM intersection file ') + wopts = C_NULL_CHAR + call shr_mpi_commrank( mpicom_CPLID, rank ) + if (rank .lt. 3) then + write(lnum,"(I0.2)")rank ! + outfile = 'intx_oa_'//trim(lnum)// '.h5m' // C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbintxoa, outfile, wopts) ! write local intx file + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing OCN-ATM intersection file ' + call shr_sys_abort(subname//' ERROR in writing OCN-ATM intersection file ') + endif endif - endif ! endif for MOABDEBUG #endif - endif ! if .not.loadmapsfromdisk - - ! we also need to compute the comm graph for the second hop, from the ocn on coupler to the - ! ocean for the intx ocean-atm context (coverage) - type1 = 3; ! fv for ocean and atm; fv-cgll does not work anyway - type2 = 3; - ierr = iMOAB_ComputeCommGraph( mboxid, mbintxoa, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & - ocn(1)%cplcompid, idintx) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for second hop, ocn-atm' - call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ocn-atm') - endif - - ! now take care of the mapper - mapper_So2a%intx_context = idintx - - if (compute_maps_online_o2a) then + ! we also need to compute the comm graph for the second hop, from the ocn on coupler to the + ! ocean for the intx ocean-atm context (coverage) + type1 = 3; ! fv for ocean and atm; fv-cgll does not work anyway + type2 = 3; + ierr = iMOAB_ComputeCommGraph( mboxid, mbintxoa, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & + ocn(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, ocn-atm' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ocn-atm') + endif if (atm_pg_active) then dm2 = "fv"//C_NULL_CHAR @@ -383,10 +369,19 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at endif else type1 = 3 ! this is type of grid, maybe should be saved on imoab app ? - call moab_map_init_rcfile( mboxid, mbaxid, mbintxoa, type1, & 'seq_maps.rc', 'ocn2atm_smapname:', 'ocn2atm_smaptype:',samegrid_ao, & wgtIdo2a, 'mapper_So2a MOAB init', esmf_map_flag) + ! need to call migrate map mesh, which will compute the cov mesh and + ! comm graph too for coverage mesh + context_id = idintx ! intx id + ierr = iMOAB_MigrateMapMesh (mboxid, mbintxoa, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, ocn(1)%cplcompid, context_id) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating ocn mesh for map ocn c2 atm ' + call shr_sys_abort(subname//' ERROR in migrating ocn mesh for map ocn c2 atm ') + endif + endif else ! samegrid_ao = TRUE @@ -543,7 +538,6 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at 'seq_maps.rc','ice2atm_smapname:','ice2atm_smaptype:',samegrid_ao, & 'mapper_Si2a initialization',esmf_map_flag, no_match) ! similar to ocn-atm mapping, do ice 2 atm mapping / set up - #ifdef HAVE_MOAB ! Call moab intx only if ATM and ICE are init in moab coupler if ((mbaxid .ge. 0) .and. (mbixid .ge. 0)) then @@ -551,7 +545,6 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at write(logunit,*) ' ' write(logunit,F00) 'Initializing MOAB mapper_Si2a' endif - appname = "ICE_ATM_COU"//C_NULL_CHAR ! idintx is a unique number of MOAB app that takes care of intx between ice and atm mesh idintx = 100*ice(1)%cplcompid + atm(1)%cplcompid ! something different, to differentiate it @@ -560,15 +553,7 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at write(logunit,*) subname,' error in registering ice atm intx' call shr_sys_abort(subname//' ERROR in registering ice atm intx') endif - - ! compute the ATM coverage mesh here on coupler pes - ! ATM mesh was redistributed to cover target (OCN) partition - ierr = iMOAB_ComputeCoverageMesh( mbixid, mbaxid, mbintxia ) - if (ierr .ne. 0) then - write(logunit,*) subname,' cannot compute source ICE coverage mesh for ATM' - call shr_sys_abort(subname//' ERROR in computing ICE-ATM coverage') - endif - + call seq_comm_getinfo(CPLID ,mpigrp=mpigrp_CPLID) if (compute_maps_online_i2a) then ierr = iMOAB_ComputeMeshIntersectionOnSphere ( mbixid, mbaxid, mbintxia ) if (ierr .ne. 0) then @@ -578,19 +563,29 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at if (iamroot_CPLID) then write(logunit,*) 'iMOAB intersection between ice and atm with id:', idintx endif - endif - - ! we also need to compute the comm graph for the second hop, from the ice on coupler to the - ! ice for the intx ice-atm context (coverage) - call seq_comm_getinfo(CPLID ,mpigrp=mpigrp_CPLID) - - type1 = 3; ! fv for ice and atm; fv-cgll does not work anyway - type2 = 3; - ierr = iMOAB_ComputeCommGraph( mbixid, mbintxia, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & - ice(1)%cplcompid, idintx) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for second hop, ice-atm' - call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ice-atm') + ! we also need to compute the comm graph for the second hop, from the ice on coupler to the + ! ice for the intx ice-atm context (coverage) +#ifdef MOABDEBUG + wopts = C_NULL_CHAR + call shr_mpi_commrank( mpicom_CPLID, rank ) + if (rank .lt. 3 ) then + write(lnum,"(I0.2)")rank ! + outfile = 'intx_ia_'//trim(lnum)// '.h5m' // C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbintxia, outfile, wopts) ! write local intx file + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing intx file ice-atm ' + call shr_sys_abort(subname//' ERROR in writing intx file ice-atm ') + endif + endif +#endif + type1 = 3; ! fv for ice and atm; fv-cgll does not work anyway + type2 = 3; + ierr = iMOAB_ComputeCommGraph( mbixid, mbintxia, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & + ice(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, ice-atm' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ice-atm') + endif endif ! now take care of the mapper mapper_Si2a%src_mbid = mbixid @@ -600,7 +595,6 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at mapper_Si2a%intx_context = idintx mapper_Si2a%weight_identifier = wgtIdi2a mapper_Si2a%mbname = 'mapper_Si2a' - if (compute_maps_online_i2a) then volumetric = 0 ! can be 1 only for FV->DGLL or FV->CGLL; if (atm_pg_active) then @@ -643,22 +637,17 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at call moab_map_init_rcfile(mbixid, mbaxid, mbintxia, type1, & 'seq_maps.rc', 'ice2atm_smapname:', 'ice2atm_smaptype:', samegrid_ao, & wgtIdi2a, 'mapper_Si2a MOAB init', esmf_map_flag) - endif - -#ifdef MOABDEBUG - wopts = C_NULL_CHAR - call shr_mpi_commrank( mpicom_CPLID, rank ) - if (rank .lt. 3 .and. compute_maps_online_i2a) then - write(lnum,"(I0.2)")rank ! - outfile = 'intx_ia_'//trim(lnum)// '.h5m' // C_NULL_CHAR - ierr = iMOAB_WriteMesh(mbintxia, outfile, wopts) ! write local intx file + context_id = idintx ! intx id + ierr = iMOAB_MigrateMapMesh (mbixid, mbintxia, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, ice(1)%cplcompid, context_id) if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing intx file ice-atm ' - call shr_sys_abort(subname//' ERROR in writing intx file ice-atm ') + write(logunit,*) subname,' error in migrating ocn mesh for map ocn c2 atm ' + call shr_sys_abort(subname//' ERROR in migrating ocn mesh for map ocn c2 atm ') endif + endif -! endif for MOABDEBUG -#endif + + endif ! if ((mbaxid .ge. 0) .and. (mbixid .ge. 0)) then ! endif for HAVE_MOAB #endif @@ -750,15 +739,7 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at mapper_Fl2a%mbname = 'mapper_Fl2a' if (.not. samegrid_al) then ! tri grid case - - ! compute the LND coverage mesh here on coupler pes - ! LND mesh was redistributed to cover target (ATM) partition - ierr = iMOAB_ComputeCoverageMesh( mblxid, mbaxid, mbintxla ) - if (ierr .ne. 0) then - write(logunit,*) subname,' cannot compute source LND coverage mesh for ATM' - call shr_sys_abort(subname//' ERROR in computing LND-ATM coverage') - endif - + call seq_comm_getinfo(CPLID ,mpigrp=mpigrp_CPLID) if (compute_maps_online_l2a) then if (iamroot_CPLID) then write(logunit,*) 'iMOAB intersection between LND and ATM with id:', idintx @@ -768,37 +749,30 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at write(logunit,*) subname,' error in computing intersection between LND and ATM' call shr_sys_abort(subname//' ERROR in computing intersection between LND and ATM') endif - endif -#ifdef MOABDEBUG - ! write intx only if true intx file: - wopts = C_NULL_CHAR - call shr_mpi_commrank( mpicom_CPLID, rank ) - if (rank .lt. 5 .and. compute_maps_online_l2a) then ! write only a few intx files - write(lnum,"(I0.2)")rank ! - outfile = 'intx_la'//trim(lnum)// '.h5m' // C_NULL_CHAR - ierr = iMOAB_WriteMesh(mbintxla, outfile, wopts) ! write local intx file + ! we also need to compute the comm graph for the second hop, from the lnd on coupler to the + ! lnd for the intx lnd-atm context (coverage) + type1 = 3; ! fv for lnd and atm; fv-cgll does not work anyway + type2 = 3; + ierr = iMOAB_ComputeCommGraph( mblxid, mbintxla, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & + lnd(1)%cplcompid, idintx) if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing intx file land atm ' - call shr_sys_abort(subname//' ERROR in writing intx file ') + write(logunit,*) subname,' error in computing comm graph for second hop, ice-atm' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ice-atm') + endif +#ifdef MOABDEBUG + ! write intx only if true intx file: + wopts = C_NULL_CHAR + call shr_mpi_commrank( mpicom_CPLID, rank ) + if (rank .lt. 5 .and. compute_maps_online_l2a) then ! write only a few intx files + write(lnum,"(I0.2)")rank ! + outfile = 'intx_la'//trim(lnum)// '.h5m' // C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbintxla, outfile, wopts) ! write local intx file + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing intx file land atm ' + call shr_sys_abort(subname//' ERROR in writing intx file ') + endif endif - endif #endif - - ! we also need to compute the comm graph for the second hop, from the lnd on coupler to the - ! lnd for the intx lnd-atm context (coverage) - call seq_comm_getinfo(CPLID ,mpigrp=mpigrp_CPLID) - - type1 = 3; ! fv for lnd and atm; fv-cgll does not work anyway - type2 = 3; - ierr = iMOAB_ComputeCommGraph( mblxid, mbintxla, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & - lnd(1)%cplcompid, idintx) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for second hop, ice-atm' - call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ice-atm') - endif - - if (compute_maps_online_l2a) then - ! need to compute weigths volumetric = 0 ! can be 1 only for FV->DGLL or FV->CGLL; if (atm_pg_active) then @@ -839,6 +813,15 @@ subroutine prep_atm_init(infodata, ocn_c2_atm, ice_c2_atm, lnd_c2_atm, iac_c2_at call moab_map_init_rcfile(mblxid, mbaxid, mbintxla, type1, & 'seq_maps.rc', 'lnd2atm_fmapname:', 'lnd2atm_fmaptype:', samegrid_al, & wgtIdl2a, 'mapper_Fl2a MOAB initialization', esmf_map_flag) + + context_id = idintx ! intx id + ierr = iMOAB_MigrateMapMesh (mblxid, mbintxla, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, lnd(1)%cplcompid, context_id) + + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating lnd mesh for map lnd c2 atm ' + call shr_sys_abort(subname//' ERROR in migrating lnd mesh for map lnd c2 atm ') + endif endif else ! the same mesh , atm and lnd use the same dofs, but restricted diff --git a/driver-moab/main/prep_lnd_mod.F90 b/driver-moab/main/prep_lnd_mod.F90 index d179456e72b9..e8a5d8faed6d 100644 --- a/driver-moab/main/prep_lnd_mod.F90 +++ b/driver-moab/main/prep_lnd_mod.F90 @@ -36,8 +36,8 @@ module prep_lnd_mod #ifdef HAVE_MOAB use iMOAB , only: iMOAB_ComputeCommGraph, iMOAB_ComputeMeshIntersectionOnSphere, & iMOAB_ComputeScalarProjectionWeights, iMOAB_DefineTagStorage, iMOAB_RegisterApplication, & - iMOAB_WriteMesh, iMOAB_GetMeshInfo, iMOAB_SetDoubleTagStorage, iMOAB_ComputeCoverageMesh, & - iMOAB_SetMapGhostLayers + iMOAB_WriteMesh, iMOAB_GetMeshInfo, iMOAB_SetDoubleTagStorage, & + iMOAB_SetMapGhostLayers, iMOAB_MigrateMapMesh use seq_comm_mct, only : num_moab_exports #endif @@ -110,6 +110,10 @@ module prep_lnd_mod ! other fields (besides frac_field and topo_field) that are mapped from glc to lnd, ! separated by elevation class character(CXX) :: glc2lnd_ec_extra_fields + +#ifdef MOABDEBUG + character*32 :: outfile, wopts, lnum +#endif !================================================================================================ contains @@ -148,7 +152,7 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln #ifdef HAVE_MOAB ! MOAB stuff integer :: ierr, idintx, rank - character*32 :: appname, outfile, wopts, lnum + character*32 :: appname character*32 :: dm1, dm2, dofnameS, dofnameT, wgtIdr2l, wgtIda2l_conservative, wgtIda2l_bilinear integer :: orderS, orderT, volumetric, noConserve, validate, fInverseDistanceMap integer :: fNoBubble, monotonicity @@ -254,12 +258,12 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln write(logunit,*) subname,' error in registering rof lnd intx' call shr_sys_abort(subname//' ERROR in registering rof lnd intx') endif + call seq_comm_getData(CPLID ,mpigrp=mpigrp_CPLID) if (samegrid_lr) then ! the same mesh , lnd and rof use the same dofs, but restricted ! we do not compute intersection, so we will have to just send data from lnd to rof and viceversa, by GLOBAL_ID matching ! so we compute just a comm graph, between lnd and rof dofs, on the coupler; target is rof ! land is full mesh - call seq_comm_getData(CPLID ,mpigrp=mpigrp_CPLID) type1 = 3; ! full mesh for lrofarnd now type2 = 3; ! fv for target land ierr = iMOAB_ComputeCommGraph( mbrxid, mblxid, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & @@ -281,15 +285,6 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln mapper_Fr2l%src_context = rof(1)%cplcompid mapper_Fr2l%intx_context = lnd(1)%cplcompid else - - ! compute the ROF coverage mesh here on coupler pes - ! ROF mesh was redistributed to cover target (LND) partition - ierr = iMOAB_ComputeCoverageMesh( mbrxid, mblxid, mbintxrl ) - if (ierr .ne. 0) then - write(logunit,*) subname,' cannot compute source ROF coverage mesh for LND' - call shr_sys_abort(subname//' ERROR in computing ROF-LND coverage') - endif - if (compute_maps_online_r2l) then ierr = iMOAB_ComputeMeshIntersectionOnSphere( mbrxid, mblxid, mbintxrl ) if (ierr .ne. 0) then @@ -300,30 +295,29 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln write(logunit,*) 'iMOAB intersection between rof and lnd with id:', idintx end if #ifdef MOABDEBUG - wopts = C_NULL_CHAR - call shr_mpi_commrank( mpicom_CPLID, rank ) - if (rank .lt. 5) then - write(lnum,"(I0.2)")rank ! - outfile = 'intx_rl_'//trim(lnum)// '.h5m' // C_NULL_CHAR - ierr = iMOAB_WriteMesh(mbintxrl, outfile, wopts) ! write local intx file - if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing intx rl file ' - call shr_sys_abort(subname//' ERROR in writing intx rl file ') - endif - endif + wopts = C_NULL_CHAR + call shr_mpi_commrank( mpicom_CPLID, rank ) + if (rank .lt. 5) then + write(lnum,"(I0.2)")rank ! + outfile = 'intx_rl_'//trim(lnum)// '.h5m' // C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbintxrl, outfile, wopts) ! write local intx file + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing intx rl file ' + call shr_sys_abort(subname//' ERROR in writing intx rl file ') + endif + endif #endif - endif + ! we also need to compute the comm graph for the second hop, from the rof on coupler to the + ! rof for the intx rof-lnd context (coverage) - ! we also need to compute the comm graph for the second hop, from the rof on coupler to the - ! rof for the intx rof-lnd context (coverage) - call seq_comm_getData(CPLID ,mpigrp=mpigrp_CPLID) - type1 = 3 ! land is FV now on coupler side - type2 = 3; - ierr = iMOAB_ComputeCommGraph( mbrxid, mbintxrl, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & - rof(1)%cplcompid, idintx) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for second hop, lnd-rof' - call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, lnd-rof') + type1 = 3 ! land is FV now on coupler side + type2 = 3; + ierr = iMOAB_ComputeCommGraph( mbrxid, mbintxrl, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & + rof(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, lnd-rof' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, lnd-rof') + endif endif ! now take care of the mapper if ( mapper_Fr2l%src_mbid .gt. -1 ) then @@ -375,10 +369,17 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln endif else type1 = 3 ! this is type of grid, maybe should be saved on imoab app ? - call moab_map_init_rcfile( mbrxid, mblxid, mbintxrl, type1, & 'seq_maps.rc', 'rof2lnd_fmapname:', 'rof2lnd_fmaptype:',samegrid_lr, & wgtIdr2l, 'mapper_Fr2l MOAB initialization', esmf_map_flag) + ! need to call migrate map mesh, which will compute the cov mesh and + ! comm graph too for coverage mesh + ierr = iMOAB_MigrateMapMesh (mbrxid, mbintxrl, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, rof(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating rof mesh for map rof c2 lnd ' + call shr_sys_abort(subname//' ERROR in migrating rof mesh for map rof c2 lnd') + endif endif endif @@ -471,27 +472,18 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln call seq_comm_getinfo(CPLID ,mpigrp=mpigrp_CPLID) if (.not. samegrid_al) then ! tri grid case - ! for bilinear maps, we need to have a layer of ghosts on source - nghlay = 1 ! number of ghost layers - nghlay_tgt = 0 - ierr = iMOAB_SetMapGhostLayers( mbintxal, nghlay, nghlay_tgt ) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in setting the number of ghost layers' - call shr_sys_abort(subname//' error in setting the number of ghost layers') - endif - - ! compute the ATM coverage mesh here on coupler pes - ! ATM mesh was redistributed to cover target (LND) partition - ierr = iMOAB_ComputeCoverageMesh( mbaxid, mblxid, mbintxal ) - if (ierr .ne. 0) then - write(logunit,*) subname,' cannot compute source ATM coverage mesh for LND' - call shr_sys_abort(subname//' ERROR in computing ATM-LND coverage') - endif - if (compute_maps_online_a2l) then if (iamroot_CPLID) then write(logunit,*) 'iMOAB intersection between atm and land with id:', idintx endif + ! for bilinear maps, we need to have a layer of ghosts on source + nghlay = 1 ! number of ghost layers + nghlay_tgt = 0 + ierr = iMOAB_SetMapGhostLayers( mbintxal, nghlay, nghlay_tgt ) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in setting the number of ghost layers' + call shr_sys_abort(subname//' error in setting the number of ghost layers') + endif ierr = iMOAB_ComputeMeshIntersectionOnSphere( mbaxid, mblxid, mbintxal ) if (ierr .ne. 0) then write(logunit,*) subname,' error in computing atm lnd intx' @@ -511,24 +503,25 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln endif endif #endif - endif - ! we also need to compute the comm graph for the second hop, from the atm on coupler to the - ! lnd for the intx atm-lnd context (coverage) - ! - if (atm_pg_active) then - type1 = 3; ! fv for atm; cgll does not work anyway - else - type1 = 1 ! this projection works (cgll to fv), but reverse does not ( fv - cgll) - endif - type2 = 3; ! land is fv in this case (separate grid) + ! we also need to compute the comm graph for the second hop, from the atm on coupler to the + ! lnd for the intx atm-lnd context (coverage) + ! + if (atm_pg_active) then + type1 = 3; ! fv for atm; cgll does not work anyway + else + type1 = 1 ! this projection works (cgll to fv), but reverse does not ( fv - cgll) + endif + type2 = 3; ! land is fv in this case (separate grid) - ierr = iMOAB_ComputeCommGraph( mbaxid, mbintxal, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & - atm(1)%cplcompid, idintx) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for second hop, atm-lnd' - call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, atm-lnd') + ierr = iMOAB_ComputeCommGraph( mbaxid, mbintxal, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & + atm(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, atm-lnd' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, atm-lnd') + endif endif + if (compute_maps_online_a2l) then volumetric = 0 ! can be 1 only for FV->DGLL or FV->CGLL; if (atm_pg_active) then @@ -584,11 +577,22 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln 'seq_maps.rc', 'atm2lnd_smapname:', 'atm2lnd_smaptype:',samegrid_al, & wgtIda2l_bilinear, 'mapper_Sa2l MOAB initialization', esmf_map_flag) + ! TODO make sure it is good enough + ! we are using the same coverage for 2 maps , one bilinear, one conservative ! read as the second one the f map, this one has the aream for land correct, so it should be fine ! the area_b for the bilinear map above is 0 ! which caused grief call moab_map_init_rcfile( mbaxid, mblxid, mbintxal, type1, & 'seq_maps.rc', 'atm2lnd_fmapname:', 'atm2lnd_fmaptype:',samegrid_al, & wgtIda2l_conservative, 'mapper_Fa2l MOAB initialization', esmf_map_flag) + ! we need to do only one map migrate, should over both maps!! + ! we have one coverage for both maps! + ierr = iMOAB_MigrateMapMesh (mbaxid, mbintxal, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, atm(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating atm mesh for map atm c2 lnd ' + call shr_sys_abort(subname//' ERROR in migrating atm mesh for map rof c2 lnd') + endif + endif else @@ -770,7 +774,6 @@ subroutine prep_lnd_mrg_moab (infodata) type(mct_aVect_sharedindices),save :: g2x_sharedindices #ifdef MOABDEBUG integer :: ierr - character*32 :: outfile, wopts, lnum #endif #ifdef MOABCOMP character(CXX) :: tagname, mct_field @@ -984,6 +987,9 @@ subroutine prep_lnd_calc_r2x_lx(timer) integer :: eri type(mct_aVect) , pointer :: r2x_rx character(*), parameter :: subname = '(prep_lnd_calc_r2x_lx)' +#ifdef MOABDEBUG + integer :: ierr +#endif !--------------------------------------------------------------- call t_drvstartf (trim(timer),barrier=mpicom_CPLID) @@ -997,6 +1003,14 @@ subroutine prep_lnd_calc_r2x_lx(timer) ! equivalent. call seq_map_map(mapper_Fr2l, r2x_rx, r2x_lx(eri), & fldlist=seq_flds_r2x_fluxes, norm=.true.) +#ifdef MOABDEBUG + if (mblxid .ge. 0 ) then ! we are on coupler pes, for sure + write(lnum,"(I0.2)")num_moab_exports + outfile = 'LndAftRofProj'//trim(lnum)//'.h5m'//C_NULL_CHAR + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR ! + ierr = iMOAB_WriteMesh(mblxid, trim(outfile), trim(wopts)) + endif +#endif enddo call t_drvstopf (trim(timer)) diff --git a/driver-moab/main/prep_ocn_mod.F90 b/driver-moab/main/prep_ocn_mod.F90 index 3f6e0f4bb291..c83b8124d22f 100644 --- a/driver-moab/main/prep_ocn_mod.F90 +++ b/driver-moab/main/prep_ocn_mod.F90 @@ -13,8 +13,6 @@ module prep_ocn_mod use seq_comm_mct, only: seq_comm_getData=>seq_comm_setptrs use seq_comm_mct, only: mboxid ! iMOAB id for mpas ocean migrated mesh to coupler pes -! use seq_comm_mct, only: mbrmapro ! iMOAB id for map read from rof2ocn map file -! use seq_comm_mct, only: mbrxoid ! iMOAB id for rof on coupler in ocean context; use seq_comm_mct, only: mbrxid ! iMOAB id of moab rof migrated to couple pes use seq_comm_mct, only: mbintxro ! iMOAB id for map read from rof2ocn map file use seq_comm_mct, only : atm_pg_active ! whether the atm uses FV mesh or not ; made true if fv_nphys > 0 @@ -209,7 +207,7 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc use iMOAB, only: iMOAB_ComputeMeshIntersectionOnSphere, iMOAB_RegisterApplication, & iMOAB_WriteMesh, iMOAB_DefineTagStorage, iMOAB_ComputeCommGraph, iMOAB_ComputeScalarProjectionWeights, & iMOAB_MigrateMapMesh, iMOAB_WriteLocalMesh, iMOAB_GetMeshInfo, iMOAB_SetDoubleTagStorage, & - iMOAB_WriteMappingWeightsToFile, iMOAB_SetMapGhostLayers, iMOAB_ComputeCoverageMesh + iMOAB_WriteMappingWeightsToFile, iMOAB_SetMapGhostLayers !--------------------------------------------------------------- ! Description ! Initialize module attribute vectors and all other non-mapping @@ -265,7 +263,7 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc integer :: rmapid, rmapid2 ! external id to identify the moab app ; 2 is for rof in ocean context (coverage) integer :: type_grid ! - integer :: context_id, direction + integer :: context_id character*32 :: prefix_output ! for writing a coverage file for debugging integer :: rank_on_cpl ! just for debugging ! these are just to zero out r2x fields on ocean @@ -439,7 +437,9 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc mapper_Fa2o%src_context = atm(1)%cplcompid mapper_Fa2o%weight_identifier = wgtIda2o_conservative mapper_Fa2o%mbname = 'mapper_Fa2o' - + ! now take care of the mapper + mapper_Fa2o%intx_context = idintx + !! updated mapper_Fa2o -- ! we also need to compute the comm graph for the second hop, from the atm on coupler to the ! atm for the intx atm-ocn context (coverage) call seq_comm_getinfo(CPLID, mpigrp=mpigrp_CPLID) @@ -447,24 +447,15 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc ! next, let us compute the ATM and OCN data transfer if (.not. samegrid_ao) then ! not a data OCN model - ! for bilinear maps, we need to have a layer of ghosts on source - nghlay = 1 ! number of ghost layers - nghlay_tgt = 0 - ierr = iMOAB_SetMapGhostLayers( mbintxao, nghlay, nghlay_tgt ) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in setting the number of ghost layers' - call shr_sys_abort(subname//' error in setting the number of ghost layers') - endif - - ! compute the ATM coverage mesh here on coupler pes - ! ATM mesh was redistributed to cover target (OCN) partition - ierr = iMOAB_ComputeCoverageMesh( mbaxid, mboxid, mbintxao ) - if (ierr .ne. 0) then - write(logunit,*) subname,' cannot compute source ATM coverage mesh for OCN' - call shr_sys_abort(subname//' ERROR in computing ATM-OCN coverage') - endif - if (compute_maps_online_a2o) then + ! for bilinear maps, we need to have a layer of ghosts on source + nghlay = 1 ! number of ghost layers + nghlay_tgt = 0 + ierr = iMOAB_SetMapGhostLayers( mbintxao, nghlay, nghlay_tgt ) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in setting the number of ghost layers' + call shr_sys_abort(subname//' error in setting the number of ghost layers') + endif ! first compute the overlap mesh between mbaxid (ATM) and mboxid (OCN) on coupler PEs ierr = iMOAB_ComputeMeshIntersectionOnSphere( mbaxid, mboxid, mbintxao ) if (ierr .ne. 0) then @@ -474,12 +465,32 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc if (iamroot_CPLID) then write(logunit,*) 'iMOAB mesh intersection completed between ATM and OCN with id:', idintx end if + if (atm_pg_active) then + type1 = 3; ! FV for ATM; CGLL does not work correctly in parallel at the moment + else + type1 = 1 ! This projection works (CGLL to FV), but reverse does not (FV - CGLL) + endif + type2 = 3; ! FV mesh on coupler OCN + ierr = iMOAB_ComputeCommGraph( mbaxid, mbintxao, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & + atm(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, ATM-OCN' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ATM-OCN') + endif +#ifdef MOABDEBUG + wopts = C_NULL_CHAR + call shr_mpi_commrank( mpicom_CPLID, rank ) + if (rank .lt. 5 .and. compute_maps_online_a2o) then + write(lnum,"(I0.2)")rank ! + outfile = 'intx_ao_'//trim(lnum)// '.h5m' // C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbintxao, outfile, wopts) ! write local intx file + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing intx file ' + call shr_sys_abort(subname//' ERROR in writing intx file ') + endif + endif +#endif end if - - ! now take care of the mapper - mapper_Fa2o%intx_context = idintx - !! updated mapper_Fa2o -- - ! To project fields from ATM to OCN grid, we need to define ! ATM a2x fields to OCN grid on coupler side tagname = trim(seq_flds_a2x_fields)//C_NULL_CHAR @@ -490,7 +501,6 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc write(logunit,*) subname,' error in defining tags for seq_flds_a2x_fields on OCN cpl' call shr_sys_abort(subname//' ERROR in coin defining tags for seq_flds_a2x_fields on OCN cpl') endif - if (compute_maps_online_a2o) then volumetric = 0 ! can be 1 only for FV->DGLL or FV->CGLL; if (atm_pg_active) then @@ -554,34 +564,18 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc call moab_map_init_rcfile(mbaxid, mboxid, mbintxao, type1, & 'seq_maps.rc', 'atm2ocn_smapname:', 'atm2ocn_smaptype:',samegrid_ao, & wgtIda2o_bilinear, 'mapper_Sa2o moab initialization', esmf_map_flag) - endif - if (atm_pg_active) then - type1 = 3; ! FV for ATM; CGLL does not work correctly in parallel at the moment - else - type1 = 1 ! This projection works (CGLL to FV), but reverse does not (FV - CGLL) - endif - type2 = 3; ! FV mesh on coupler OCN - ierr = iMOAB_ComputeCommGraph( mbaxid, mbintxao, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & - atm(1)%cplcompid, idintx) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for second hop, ATM-OCN' - call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ATM-OCN') - endif + context_id = idintx + ! again, one coverage set and coverage graph for 2 different maps + ierr = iMOAB_MigrateMapMesh (mbaxid, mbintxao, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, atm(1)%cplcompid, context_id) -#ifdef MOABDEBUG - wopts = C_NULL_CHAR - call shr_mpi_commrank( mpicom_CPLID, rank ) - if (rank .lt. 5 .and. compute_maps_online_a2o) then - write(lnum,"(I0.2)")rank ! - outfile = 'intx_ao_'//trim(lnum)// '.h5m' // C_NULL_CHAR - ierr = iMOAB_WriteMesh(mbintxao, outfile, wopts) ! write local intx file if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing intx file ' - call shr_sys_abort(subname//' ERROR in writing intx file ') + write(logunit,*) subname,' error in migrating atm mesh for map atm c2 ocn ' + call shr_sys_abort(subname//' ERROR in migrating atm mesh for map atm c2 ocn ') endif + endif -#endif else ! if (samegrid_ao) ! ATM and OCN components use the same mesh and DoF numbering (OCN is a subset of ATM); @@ -748,27 +742,6 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc write(logunit,*) subname,' error in registering rof 2 ocn moab map ' call shr_sys_abort(subname//' ERROR in registering rof 2 ocn moab map ') endif - ! integer, public :: mboxid ! iMOAB id for mpas ocean already migrated mesh to coupler pes - ! this is a special rof mesh redistribution, for the ocean context - ! it will be used to project from rof to ocean - ! the mesh will be migrated, to be able to do the second hop - ! appname = "ROF_OCOU"//C_NULL_CHAR - ! ! rmapid is a unique external number of MOAB app that identifies runoff on coupler side - ! rmapid2 = 100*rof(1)%cplcompid ! this is a special case, because we also have a regular coupler instance mbrxid - ! ierr = iMOAB_RegisterApplication(trim(appname), mpicom_CPLID, rmapid2, mbrxoid) - ! if (ierr .ne. 0) then - ! write(logunit,*) subname,' error in registering rof on coupler in ocean context ' - ! call shr_sys_abort(subname//' ERROR in registering rof on coupler in ocean context ') - ! endif - ! this code was moved from prep_rof_ocn_moab, because we will do everything on coupler side, not - ! needed to be on joint comm anymore for the second hop - - ! need to compute coverage of rof over ocean, and comm graph for sending from rof to rof over ocean - ierr = iMOAB_ComputeCoverageMesh( mbrxid, mboxid, mbintxro ) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in compute coverage mesh rof for ocean' - call shr_sys_abort(subname//' ERROR in compute coverage mesh rof for ocean ') - endif if (iamroot_CPLID) then write(logunit,*) ' ' write(logunit,F00) 'Initializing MOAB mapper_Rr2o_liq' @@ -779,44 +752,20 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc 'seq_maps.rc', 'rof2ocn_liq_rmapname:', 'rof2ocn_liq_rmaptype:',samegrid_ro, & wgtIdr2o_conservative, 'mapper_Rr2o_liq moab initialization',esmf_map_flag) - ierr = iMOAB_ComputeCommGraph( mbrxid, mbintxro, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type_grid, & - type_grid, rof(1)%cplcompid, rmapid ) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in compute graph rof - rof cov for ocean' - call shr_sys_abort(subname//' ERROR in compute graph rof - rof cov for ocean ') - endif - ! it read on the coupler side, from file, the scrip mosart, that has a full mesh; - ! also migrate rof mesh on coupler pes, in ocean context, mbrxoid (this will be like coverage mesh, - ! it will cover ocean target per process) - ! map between rof 2 ocn is in mbrmapro ; - ! after this, the sending of tags for second hop (ocn context) will use the new par comm graph, - ! that has more precise info, that got created - - ! type1 = 3 ! fv mesh nowadays - ! direction = 1 ! - ! context_id = ocn(1)%cplcompid - ! this creates a par comm graph between mbrxid and mbrxoid, with ids rof(1)%cplcompid, context ocn(1)%cplcompid + type1 = 3 ! fv mesh nowadays + context_id = rmapid ! ocn(1)%cplcompid + ! this creates a par comm graph between mbrxid and mbintxro, with ids rof(1)%cplcompid, rmapid (rmapid is 100*src+tgt) ! this will be used in send/receive mappers - ! ierr = iMOAB_MigrateMapMesh (mbrxid, mbrmapro, mbrxoid, mpicom_CPLID, mpigrp_CPLID, & - ! mpigrp_CPLID, type1, rof(1)%cplcompid, context_id, direction) + ierr = iMOAB_MigrateMapMesh (mbrxid, mbintxro, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, rof(1)%cplcompid, context_id) - ! if (ierr .ne. 0) then - ! write(logunit,*) subname,' error in migrating rof mesh for map rof c2 ocn ' - ! call shr_sys_abort(subname//' ERROR in migrating rof mesh for map rof c2 ocn ') - ! endif + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating rof mesh for map rof c2 ocn ' + call shr_sys_abort(subname//' ERROR in migrating rof mesh for map rof c2 ocn ') + endif ! if (iamroot_CPLID) then ! write(logunit,*) subname,' migrated mesh for map rof 2 ocn ' ! endif - if (mbrxid .ge. 0) then ! we are on coupler side pes - tagname=trim(seq_flds_r2x_fields)//C_NULL_CHAR - tagtype = 1 ! dense, double - numco = 1 ! only 1 component DoF per node - ierr = iMOAB_DefineTagStorage(mbrxid, tagname, tagtype, numco, tagindex ) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in defining ' // trim(seq_flds_r2x_fields) // ' tags on coupler side in MOAB' - call shr_sys_abort(subname//' ERROR in defining MOAB tags ') - endif - endif if (mboxid .ge. 0) then ! we are on coupler side pes, for ocean mesh tagname=trim(seq_flds_r2x_fields)//C_NULL_CHAR @@ -854,19 +803,6 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc endif deallocate (tmparray) - ! now we have to populate the map with the right moab attributes, so that it does the right projection -#ifdef MOABDEBUG - if (mbrxid.ge.0) then ! we are on coupler PEs - call mpi_comm_rank(mpicom_CPLID, rank_on_cpl , ierr) - if (rank_on_cpl .lt. 4) then - prefix_output = "rof_cov"//CHAR(0) - ierr = iMOAB_WriteLocalMesh(mbrxid, prefix_output) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing coverage mesh rof 2 ocn ' - endif - endif - endif -#endif ! now take care of the mapper for MOAB mapper_Rr2o_liq if ( mapper_Rr2o_liq%src_mbid .gt. -1 ) then if (iamroot_CPLID) then @@ -875,7 +811,7 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc endif endif mapper_Rr2o_liq%src_mbid = mbrxid - mapper_Rr2o_liq%tgt_mbid = mboxid ! this is special, it will really need this coverage type mesh + mapper_Rr2o_liq%tgt_mbid = mboxid ! this is similar to a regular intx scenario mapper_Rr2o_liq%intx_mbid = mbintxro mapper_Rr2o_liq%src_context = rof(1)%cplcompid !mapper_Rr2o_liq%intx_context = ocn(1)%cplcompid ! this context was used in migrate mesh @@ -902,14 +838,12 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc write(logunit,F00) 'Initializing MOAB mapper_Rr2o_ice same as mapper_Rr2o_liq' end if mapper_Rr2o_ice%src_mbid = mbrxid - mapper_Rr2o_ice%tgt_mbid = mboxid ! special + mapper_Rr2o_ice%tgt_mbid = mboxid mapper_Rr2o_ice%intx_mbid = mbintxro mapper_Rr2o_ice%src_context = rof(1)%cplcompid - !mapper_Rr2o_ice%intx_context = ocn(1)%cplcompid ! this context was used in migrate mesh mapper_Rr2o_ice%intx_context = rmapid ! read map is the same context as intersection now mapper_Rr2o_ice%weight_identifier = wgtIdr2o_conservative mapper_Rr2o_ice%mbname = 'mapper_Rr2o_ice' - ! mapper_Rr2o_ice%read_map = .true. #endif if (flood_present) then if (iamroot_CPLID) then @@ -926,10 +860,9 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc write(logunit,F00) 'Initializing MOAB mapper_Fr2o' end if mapper_Fr2o%src_mbid = mbrxid - mapper_Fr2o%tgt_mbid = mboxid ! special + mapper_Fr2o%tgt_mbid = mboxid mapper_Fr2o%intx_mbid = mbintxro mapper_Fr2o%src_context = rof(1)%cplcompid - !mapper_Fr2o%intx_context = ocn(1)%cplcompid ! this context was used in migrate mesh mapper_Fr2o%intx_context = rmapid ! read map is the same context as intersection now mapper_Fr2o%weight_identifier = wgtIdr2o_conservative mapper_Fr2o%mbname = 'mapper_Fr2o' @@ -2943,7 +2876,14 @@ subroutine prep_ocn_calc_r2x_ox(timer) type(mct_avect), pointer :: r2x_rx character(*), parameter :: subname = '(prep_ocn_calc_r2x_ox)' !--------------------------------------------------------------- - +#ifdef MOABDEBUG + if (mboxid .ge. 0 ) then ! we are on coupler pes, for sure + write(lnum,"(I0.2)")num_moab_exports + outfile = 'OcnCpl_Bef_r2x_ox_'//trim(lnum)//'.h5m'//C_NULL_CHAR + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(mboxid, trim(outfile), trim(wopts)) + endif +#endif call t_drvstartf (trim(timer),barrier=mpicom_CPLID) do eri = 1,num_inst_rof r2x_rx => component_get_c2x_cx(rof(eri)) @@ -2957,6 +2897,14 @@ subroutine prep_ocn_calc_r2x_ox(timer) endif enddo call t_drvstopf (trim(timer)) +#ifdef MOABDEBUG + if (mboxid .ge. 0 ) then ! we are on coupler pes, for sure + write(lnum,"(I0.2)")num_moab_exports + outfile = 'OcnCpl_r2x_ox_'//trim(lnum)//'.h5m'//C_NULL_CHAR + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(mboxid, trim(outfile), trim(wopts)) + endif +#endif end subroutine prep_ocn_calc_r2x_ox diff --git a/driver-moab/main/prep_rof_mod.F90 b/driver-moab/main/prep_rof_mod.F90 index 007ff419d7a4..3dc96cd294a6 100644 --- a/driver-moab/main/prep_rof_mod.F90 +++ b/driver-moab/main/prep_rof_mod.F90 @@ -166,7 +166,7 @@ subroutine prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) use iMOAB, only: iMOAB_ComputeMeshIntersectionOnSphere, iMOAB_RegisterApplication, & iMOAB_WriteMesh, iMOAB_DefineTagStorage, iMOAB_ComputeCommGraph, & - iMOAB_ComputeScalarProjectionWeights, iMOAB_GetMeshInfo, iMOAB_ComputeCoverageMesh + iMOAB_ComputeScalarProjectionWeights, iMOAB_GetMeshInfo, iMOAB_MigrateMapMesh !--------------------------------------------------------------- ! Description ! Initialize module attribute vectors and all other non-mapping @@ -372,15 +372,6 @@ subroutine prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) mapper_Fl2r%src_context = lnd(1)%cplcompid mapper_Fl2r%intx_context = rof(1)%cplcompid else - - ! compute the LND coverage mesh here on coupler pes - ! LND mesh was redistributed to cover target (ROF) partition - ierr = iMOAB_ComputeCoverageMesh( mblxid, mbrxid, mbintxlr ) - if (ierr .ne. 0) then - write(logunit,*) subname,' cannot compute source LND coverage mesh for ROF' - call shr_sys_abort(subname//' ERROR in computing LND-ROF coverage') - endif - ! if we are not loading maps from disk, compute the intersection mesh between ! LND and ROF meshes if (compute_maps_online_l2r) then @@ -393,30 +384,29 @@ subroutine prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) write(logunit,*) 'iMOAB intersection between LND and ROF with id:', idintx end if #ifdef MOABDEBUG - wopts = C_NULL_CHAR - call shr_mpi_commrank( mpicom_CPLID, rank ) - if (rank .lt. 3 .and. compute_maps_online_l2r) then - write(lnum,"(I0.2)")rank ! - outfile = 'intx_lr_'//trim(lnum)// '.h5m' // C_NULL_CHAR - ierr = iMOAB_WriteMesh(mbintxlr, outfile, wopts) ! write local intx file - if (ierr .ne. 0) then - write(logunit,*) subname,' error in writing LND-ROF intersection mesh file ' - call shr_sys_abort(subname//' ERROR in writing LND-ROF intersection mesh file ') + wopts = C_NULL_CHAR + call shr_mpi_commrank( mpicom_CPLID, rank ) + if (rank .lt. 3 .and. compute_maps_online_l2r) then + write(lnum,"(I0.2)")rank ! + outfile = 'intx_lr_'//trim(lnum)// '.h5m' // C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbintxlr, outfile, wopts) ! write local intx file + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing LND-ROF intersection mesh file ' + call shr_sys_abort(subname//' ERROR in writing LND-ROF intersection mesh file ') + endif endif - endif #endif + ! we also need to compute the comm graph for the second hop, from the lnd on coupler to the + ! lnd for the intx LND-ROF context (coverage) + type1 = 3 ! land is FV now on coupler side + type2 = 3; + ierr = iMOAB_ComputeCommGraph( mblxid, mbintxlr, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & + lnd(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, LND-ROF' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, LND-ROF') + endif end if - - ! we also need to compute the comm graph for the second hop, from the lnd on coupler to the - ! lnd for the intx LND-ROF context (coverage) - type1 = 3 ! land is FV now on coupler side - type2 = 3; - ierr = iMOAB_ComputeCommGraph( mblxid, mbintxlr, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & - lnd(1)%cplcompid, idintx) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for second hop, LND-ROF' - call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, LND-ROF') - endif ! now take care of the mapper if ( mapper_Fl2r%src_mbid .gt. -1 ) then if (iamroot_CPLID) then @@ -469,6 +459,15 @@ subroutine prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) call moab_map_init_rcfile( mblxid, mbrxid, mbintxlr, type1, & 'seq_maps.rc', 'lnd2rof_fmapname:', 'lnd2rof_fmaptype:',samegrid_lr, & wgtIdl2r, 'mapper_Fl2r MOAB initialization', esmf_map_flag) + + ! this creates a par comm graph between mblxid and mbintxlr, with ids lnd(1)%cplcompid + ierr = iMOAB_MigrateMapMesh (mblxid, mbintxlr, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, lnd(1)%cplcompid, idintx) + + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating lnd mesh for map lnd c2 rof ' + call shr_sys_abort(subname//' ERROR in migrating lnd mesh for map lnd c2 rof ') + endif end if end if ! if ((mblxid .ge. 0) .and. (mbrxid .ge. 0)) endif ! samegrid_lr @@ -563,14 +562,6 @@ subroutine prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) call shr_sys_abort(subname//' ERROR in registering ATM-ROF mesh intersection context') endif - ! compute the OCN coverage mesh here on coupler pes - ! ATM mesh was redistributed to cover target (ROF) partition - ierr = iMOAB_ComputeCoverageMesh( mbaxid, mbrxid, mbintxar ) - if (ierr .ne. 0) then - write(logunit,*) subname,' cannot compute source ATM coverage mesh for ROF' - call shr_sys_abort(subname//' ERROR in computing ATM-ROF coverage') - endif - ! if we are not loading maps from disk, compute the intersection mesh between ! ATM and ROF meshes if (compute_maps_online_a2r) then @@ -595,23 +586,22 @@ subroutine prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) endif endif #endif + ! we also need to compute the comm graph for the second hop, from the atm on coupler to the + ! atm for the intersection of ATM-ROF context (coverage) + call seq_comm_getData(CPLID ,mpigrp=mpigrp_CPLID) + if (atm_pg_active) then + type1 = 3; ! FV for both rof and atm; FV-CGLL does not work anyway + else + type1 = 1 ! this does not work anyway in this direction + endif + type2 = 3; + ierr = iMOAB_ComputeCommGraph( mbaxid, mbintxar, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & + atm(1)%cplcompid, idintx) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in computing comm graph for second hop, ATM-ROF' + call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ATM-ROF') + endif end if - - ! we also need to compute the comm graph for the second hop, from the atm on coupler to the - ! atm for the intersection of ATM-ROF context (coverage) - call seq_comm_getData(CPLID ,mpigrp=mpigrp_CPLID) - if (atm_pg_active) then - type1 = 3; ! FV for both rof and atm; FV-CGLL does not work anyway - else - type1 = 1 ! this does not work anyway in this direction - endif - type2 = 3; - ierr = iMOAB_ComputeCommGraph( mbaxid, mbintxar, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, type1, type2, & - atm(1)%cplcompid, idintx) - if (ierr .ne. 0) then - write(logunit,*) subname,' error in computing comm graph for second hop, ATM-ROF' - call shr_sys_abort(subname//' ERROR in computing comm graph for second hop, ATM-ROF') - endif ! now take care of the mapper if ( mapper_Fa2r%src_mbid .gt. -1 ) then if (iamroot_CPLID) then @@ -677,10 +667,17 @@ subroutine prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) else type1 = 3 ! this is type of grid, maybe should be saved on imoab app ? - call moab_map_init_rcfile( mbaxid, mbrxid, mbintxar, type1, & 'seq_maps.rc', 'atm2rof_fmapname:', 'atm2rof_fmaptype:',samegrid_ar, & wgtIda2r, 'mapper_Fa2r MOAB initialization', esmf_map_flag) + ! this creates a par comm graph between mblxid and mbintxlr, with ids lnd(1)%cplcompid + ierr = iMOAB_MigrateMapMesh (mbaxid, mbintxar, mpicom_CPLID, mpigrp_CPLID, & + mpigrp_CPLID, type1, atm(1)%cplcompid, idintx) + + if (ierr .ne. 0) then + write(logunit,*) subname,' error in migrating atm mesh for map atm c2 rof ' + call shr_sys_abort(subname//' ERROR in migrating atm mesh for map atm c2 rof ') + endif end if end if ! if ((mbrxid .ge. 0) .and. (mbaxid .ge. 0)) diff --git a/driver-moab/main/seq_flux_mct.F90 b/driver-moab/main/seq_flux_mct.F90 index 018d128b5376..f5f7ff1b07e1 100644 --- a/driver-moab/main/seq_flux_mct.F90 +++ b/driver-moab/main/seq_flux_mct.F90 @@ -191,6 +191,7 @@ module seq_flux_mct integer :: index_xao_So_ssq integer :: index_xao_So_duu10n integer :: index_xao_So_u10 + integer :: index_xao_So_u10withgusts integer :: index_xao_So_fswpen integer :: index_xao_So_warm_diurn integer :: index_xao_So_salt_diurn @@ -1510,6 +1511,7 @@ subroutine seq_flux_atmocn_mct(infodata, tod, dt, a2x, o2x, xao) index_xao_So_re = mct_aVect_indexRA(xao,'So_re') index_xao_So_ssq = mct_aVect_indexRA(xao,'So_ssq') index_xao_So_u10 = mct_aVect_indexRA(xao,'So_u10') + index_xao_So_u10withgusts = mct_aVect_indexRA(xao,'So_u10withgusts') index_xao_So_duu10n = mct_aVect_indexRA(xao,'So_duu10n') index_xao_Faox_taux = mct_aVect_indexRA(xao,'Faox_taux') index_xao_Faox_tauy = mct_aVect_indexRA(xao,'Faox_tauy') @@ -1787,6 +1789,7 @@ subroutine seq_flux_atmocn_mct(infodata, tod, dt, a2x, o2x, xao) xao%rAttr(index_xao_Faox_lwup,n) = lwup(n) xao%rAttr(index_xao_So_duu10n,n) = duu10n(n) xao%rAttr(index_xao_So_u10 ,n) = sqrt(duu10n(n)) + xao%rAttr(index_xao_So_u10withgusts,n) = sqrt(duu10n(n)) xao%rAttr(index_xao_So_warm_diurn ,n) = warm(n) xao%rAttr(index_xao_So_salt_diurn ,n) = salt(n) xao%rAttr(index_xao_So_speed_diurn ,n) = speed(n) diff --git a/driver-moab/main/seq_frac_mct.F90 b/driver-moab/main/seq_frac_mct.F90 index 031845a102ef..ba892ce63c33 100644 --- a/driver-moab/main/seq_frac_mct.F90 +++ b/driver-moab/main/seq_frac_mct.F90 @@ -140,7 +140,7 @@ module seq_frac_mct ! !USES: - use shr_kind_mod , only: R8 => SHR_KIND_R8 + use shr_kind_mod , only: R8 => SHR_KIND_R8, IN=>SHR_KIND_IN use shr_sys_mod use shr_const_mod @@ -287,6 +287,10 @@ subroutine seq_frac_init( infodata, & !----- local ----- type(mct_ggrid), pointer :: dom_a + type(mct_gsmap), pointer :: gsmap_a ! see if we can get from here the global ids (missing from dom_a) + type(mct_gsmap), pointer :: gsmap_l + type(mct_gsmap), pointer :: gsmap_r + type(mct_gsmap), pointer :: gsmap_i ! ofrac on ice error type(mct_ggrid), pointer :: dom_i type(mct_ggrid), pointer :: dom_l type(mct_ggrid), pointer :: dom_o @@ -326,11 +330,13 @@ subroutine seq_frac_init( infodata, & character*32 :: wgtIdef real(r8), allocatable :: tagValues(:) ! used for setting some default tags integer , allocatable :: GlobalIds(:) ! used for setting values associated with ids + integer(IN), pointer :: dof(:) ! integer nvert(3), nvise(3), nbl(3), nsurf(3), nvisBC(3) integer kgg ! index in global number attribute, used for global id in MOAB integer idintx ! used for context for intx atm - ocn integer id_join ! used for example for atm%cplcompid integer :: mpicom ! we are on coupler PES here + integer :: my_task ! character(30) :: outfile, wopts @@ -361,6 +367,7 @@ subroutine seq_frac_init( infodata, & dom_w => component_get_dom_cx(wav) dom_z => component_get_dom_cx(iac) + mpicom = seq_comm_mpicom(CPLID) debug_old = seq_frac_debug seq_frac_debug = 2 @@ -469,17 +476,21 @@ subroutine seq_frac_init( infodata, & tagname = 'lfrin'//C_NULL_CHAR ! 'lfrin' allocate(tagValues(lSize) ) tagValues = dom_l%data%rAttr(kf,:) - kgg = mct_aVect_indexIA(dom_l%data ,"GlobGridNum" ,perrWith=subName) - allocate(GlobalIds(lSize)) - GlobalIds = dom_l%data%iAttr(kgg,:) - + !kgg = mct_aVect_indexIA(dom_l%data ,"GlobGridNum" ,perrWith=subName) + !allocate(GlobalIds(lSize)) + !GlobalIds = dom_l%data%iAttr(kgg,:) + gsmap_l => component_get_gsmap_cx(lnd) ! gsmap_lx + call mpi_comm_rank(mpicom,my_task,ierr) + ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid + call mct_gsMap_orderedPoints(gsmap_l, my_task, dof) ! ent_type should be 3, FV - ierr = iMOAB_SetDoubleTagStorageWithGid ( mblxid, tagname, lSize , ent_type, tagValues, GlobalIds ) + ierr = iMOAB_SetDoubleTagStorageWithGid ( mblxid, tagname, lSize , ent_type, tagValues, dof ) if (ierr .ne. 0) then write(logunit,*) subname,' error in setting lfrin on lnd ' call shr_sys_abort(subname//' ERROR in setting lfrin on lnd') endif - deallocate(GlobalIds) + !deallocate(GlobalIds) + deallocate(dof) deallocate(tagValues) endif @@ -526,16 +537,21 @@ subroutine seq_frac_init( infodata, & tagname = 'rfrac'//C_NULL_CHAR ! 'rfrac' allocate(tagValues(lSize) ) tagValues = dom_r%data%rAttr(kf,:) - kgg = mct_aVect_indexIA(dom_r%data ,"GlobGridNum" ,perrWith=subName) - allocate(GlobalIds(lSize)) - GlobalIds = dom_r%data%iAttr(kgg,:) + !kgg = mct_aVect_indexIA(dom_r%data ,"GlobGridNum" ,perrWith=subName) + !allocate(GlobalIds(lSize)) + !GlobalIds = dom_r%data%iAttr(kgg,:) + gsmap_r => component_get_gsmap_cx(rof) ! gsmap_rx + call mpi_comm_rank(mpicom,my_task,ierr) + ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid + call mct_gsMap_orderedPoints(gsmap_r, my_task, dof) ! again, we are setting on the river instance that is also used for ocean coupling - ierr = iMOAB_SetDoubleTagStorageWithGid ( mbrxid, tagname, lSize , ent_type, tagValues, GlobalIds ) + ierr = iMOAB_SetDoubleTagStorageWithGid ( mbrxid, tagname, lSize , ent_type, tagValues, dof ) if (ierr .ne. 0) then write(logunit,*) subname,' error in setting rfrac on rof ' call shr_sys_abort(subname//' ERROR in setting rfrac on rof ') endif - deallocate(GlobalIds) + !deallocate(GlobalIds) + deallocate(dof) deallocate(tagValues) endif @@ -590,23 +606,49 @@ subroutine seq_frac_init( infodata, & tagname = 'ofrac'//C_NULL_CHAR ! 'ofrac' allocate(tagValues(lSize) ) tagValues = dom_i%data%rAttr(kf,:) - kgg = mct_aVect_indexIA(dom_i%data ,"GlobGridNum" ,perrWith=subName) - allocate(GlobalIds(lSize)) - GlobalIds = dom_i%data%iAttr(kgg,:) - - ierr = iMOAB_SetDoubleTagStorageWithGid ( mbixid, tagname, lSize , ent_type, tagValues, GlobalIds ) + !kgg = mct_aVect_indexIA(dom_i%data ,"GlobGridNum" ,perrWith=subName) + !allocate(GlobalIds(lSize)) + !GlobalIds = dom_i%data%iAttr(kgg,:) + gsmap_i => component_get_gsmap_cx(ice) ! gsmap_ix + call mpi_comm_rank(mpicom,my_task,ierr) + ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid + call mct_gsMap_orderedPoints(gsmap_i, my_task, dof) + ierr = iMOAB_SetDoubleTagStorageWithGid ( mbixid, tagname, lSize , ent_type, tagValues, dof ) if (ierr .ne. 0) then write(logunit,*) subname,' error in setting ofrac on ice ' call shr_sys_abort(subname//' ERROR in setting ofrac on ice ') endif - deallocate(GlobalIds) + !deallocate(GlobalIds) deallocate(tagValues) + deallocate(dof) +#ifdef MOABDEBUG + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR + if (mbixid .ge. 0 ) then + outfile = 'iceCplInit1Fr.h5m'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbixid, trim(outfile), trim(wopts)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing mesh ' + call shr_sys_abort(subname//' ERROR in writing mesh ') + endif + endif +#endif endif if (atm_present) then mapper_i2a => prep_atm_get_mapper_Fi2a() call seq_map_map(mapper_i2a,fractions_i,fractions_a,fldlist='ofrac',norm=.false.) endif +#ifdef MOABDEBUG + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR + if (mbaxid .ge. 0 ) then + outfile = 'atmCplInit1Fr.h5m'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbaxid, trim(outfile), trim(wopts)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing mesh ' + call shr_sys_abort(subname//' ERROR in writing mesh ') + endif + endif +#endif end if ! end of ice_present @@ -674,6 +716,17 @@ subroutine seq_frac_init( infodata, & deallocate(tagValues) mapper_o2a => prep_atm_get_mapper_Fo2a() call seq_map_map(mapper_o2a, fractions_o, fractions_a, fldlist='ofrac',norm=.false.) +#ifdef MOABDEBUG + wopts = ';PARALLEL=WRITE_PART'//C_NULL_CHAR + if (mbaxid .ge. 0 ) then + outfile = 'atmCplInit2Fr.h5m'//C_NULL_CHAR + ierr = iMOAB_WriteMesh(mbaxid, trim(outfile), trim(wopts)) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in writing mesh ' + call shr_sys_abort(subname//' ERROR in writing mesh ') + endif + endif +#endif endif if (atm_present) then @@ -715,16 +768,21 @@ subroutine seq_frac_init( infodata, & tagname = 'lfrac'//C_NULL_CHAR ! 'lfrac allocate(tagValues(lSize) ) tagValues = fractions_a%rAttr(kl,:) - kgg = mct_aVect_indexIA(dom_a%data ,"GlobGridNum" ,perrWith=subName) - allocate(GlobalIds(lSize)) - GlobalIds = dom_a%data%iAttr(kgg,:) + !kgg = mct_aVect_indexIA(dom_a%data ,"GlobGridNum" ,perrWith=subName) + !allocate(GlobalIds(lSize)) + !GlobalIds = dom_a%data%iAttr(kgg,:) ! set on atmosphere instance - ierr = iMOAB_SetDoubleTagStorageWithGid ( mbaxid, tagname, lSize , ent_type, tagValues, GlobalIds ) + gsmap_a => component_get_gsmap_cx(atm) ! gsmap_ax + call mpi_comm_rank(mpicom,my_task,ierr) + ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid + call mct_gsMap_orderedPoints(gsmap_a, my_task, dof) + ierr = iMOAB_SetDoubleTagStorageWithGid ( mbaxid, tagname, lSize , ent_type, tagValues, dof ) if (ierr .ne. 0) then write(logunit,*) subname,' error in setting lfrac on atm ' call shr_sys_abort(subname//' ERROR in setting lfrac on atm ') endif - deallocate(GlobalIds) + !deallocate(GlobalIds) + deallocate(dof) deallocate(tagValues) endif else if (lnd_present) then @@ -737,6 +795,38 @@ subroutine seq_frac_init( infodata, & if (atm_frac_correct) fractions_a%rAttr(kl,n) = 1.0_r8 endif enddo + ! case --res r05_r05 --compset RMOSGPCC is coming here + ! there are no active ice or ocn comps + ! ist is almost like above, except the ofrac is modified too + if (mbaxid .ge. 0 ) then ! // + tagname = 'lfrac'//C_NULL_CHAR ! 'lfrac + allocate(tagValues(lSize) ) + tagValues = fractions_a%rAttr(kl,:) + kgg = mct_aVect_indexIA(dom_a%data ,"GlobGridNum" ,perrWith=subName) + !allocate(GlobalIds(lSize)) + ! GlobalIds = dom_a%data%iAttr(kgg,:) + gsmap_a => component_get_gsmap_cx(atm) ! gsmap_ax + call mpi_comm_rank(mpicom,my_task,ierr) + ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid + call mct_gsMap_orderedPoints(gsmap_a, my_task, dof) + ! set on atmosphere instance + ierr = iMOAB_SetDoubleTagStorageWithGid ( mbaxid, tagname, lSize , ent_type, tagValues, dof ) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in setting lfrac on atm ' + call shr_sys_abort(subname//' ERROR in setting lfrac on atm ') + endif + ! unlike above, set ofrac too + tagname = 'ofrac'//C_NULL_CHAR ! 'ofrac + tagValues = fractions_a%rAttr(ko,:) + ! set on atmosphere instance, ofrac value + ierr = iMOAB_SetDoubleTagStorageWithGid ( mbaxid, tagname, lSize , ent_type, tagValues, dof ) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in setting ofrac on atm ' + call shr_sys_abort(subname//' ERROR in setting ofrac on atm ') + endif + deallocate(dof) + deallocate(tagValues) + endif endif endif @@ -753,6 +843,24 @@ subroutine seq_frac_init( infodata, & kk = mct_aVect_indexRA(fractions_l,"lfrin",perrWith=subName) kl = mct_aVect_indexRA(fractions_l,"lfrac",perrWith=subName) fractions_l%rAttr(kl,:) = fractions_l%rAttr(kk,:) + ! still need to do MOAB case, basically copy the value of lfrin to lfrac, on land + ! fractions / on land moab instance; otherwise would stay at 0 :) + ierr = iMOAB_GetMeshInfo ( mblxid, nvert, nvise, nbl, nsurf, nvisBC ) + arrSize = nvise(1) ! there is one tag that we need to copy + allocate(tagValues(arrSize) ) + tagname = 'lfrin'//C_NULL_CHAR + ierr = iMOAB_GetDoubleTagStorage ( mbixid, tagname, arrSize , ent_type, tagValues) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in getting lfrin on lnd ' + call shr_sys_abort(subname//' ERROR in getting lfrin on lnd ') + endif + tagname = 'lfrac'//C_NULL_CHAR + ierr = iMOAB_SetDoubleTagStorage ( mbixid, tagname, arrSize , ent_type, tagValues) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in setting lfrac on lnd ' + call shr_sys_abort(subname//' ERROR in setting lfrac on lnd ') + endif + deallocate(tagValues) end if end if if (lnd_present .and. rof_present) then @@ -849,22 +957,25 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ !----- local ----- type(mct_aVect), pointer :: i2x_i type(mct_ggrid), pointer :: dom_i - type(mct_ggrid), pointer :: dom_o ! introduced just to update ocean moab fractions + !type(mct_ggrid), pointer :: dom_o ! introduced just to update ocean moab fractions + type(mct_gsmap), pointer :: gsmap_i ! ofrac on ice error logical :: atm_present ! true => atm is present logical :: ice_present ! true => ice is present logical :: ocn_present ! true => ocn is present integer :: n integer :: ki, kl, ko, kf real(r8),allocatable :: fcorr(:) + integer(IN), pointer, save :: dof(:) ! logical, save :: first_time = .true. ! moab integer :: ierr, kgg integer , save :: lSize, ent_type + integer mpicom, my_task character(CXX) :: tagname real(r8), allocatable, save :: tagValues(:) ! used for setting some tags real(r8), allocatable, save :: tagValuesOfrac(:) ! used for setting some tags - integer , allocatable, save :: GlobalIds(:) ! used for setting values associated with ids + !integer , allocatable, save :: GlobalIds(:) ! used for setting values associated with ids #ifdef MOABDEBUG character(len=100) :: outfile, wopts, lnum #endif @@ -891,8 +1002,9 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ dom_i => component_get_dom_cx(ice) i2x_i => component_get_c2x_cx(ice) + mpicom = seq_comm_mpicom(CPLID) - dom_o => component_get_dom_cx(ocn) ! + !dom_o => component_get_dom_cx(ocn) ! if (ice_present) then call mct_aVect_copy(i2x_i, fractions_i, "Si_ifrac", "ifrac") @@ -908,9 +1020,13 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ if (first_time) then ! allocate some local arrays lSize = mct_aVect_lSize(dom_i%data) allocate(tagValues(lSize) ) - allocate(GlobalIds(lSize) ) - kgg = mct_aVect_indexIA(dom_o%data ,"GlobGridNum" ,perrWith=subName) - GlobalIds = dom_i%data%iAttr(kgg,:) + !allocate(GlobalIds(lSize) ) + !kgg = mct_aVect_indexIA(dom_o%data ,"GlobGridNum" ,perrWith=subName) + !GlobalIds = dom_i%data%iAttr(kgg,:) + gsmap_i => component_get_gsmap_cx(ice) ! gsmap_ix + call mpi_comm_rank(mpicom,my_task,ierr) + ! Determine global gridpoint number attribute, GlobGridNum, automatically in ggrid + call mct_gsMap_orderedPoints(gsmap_i, my_task, dof) allocate (tagValuesOfrac(local_size_mb_ocn)) ent_type = 1 ! cells for mpas sea ice first_time = .false. @@ -923,7 +1039,7 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ ! fraclist_i = 'afrac:ifrac:ofrac' ! tagValues = fractions_i%rAttr(ki,:) - ierr = iMOAB_SetDoubleTagStorageWithGid ( mbixid, tagname, lSize , ent_type, tagValues, GlobalIds ) + ierr = iMOAB_SetDoubleTagStorageWithGid ( mbixid, tagname, lSize , ent_type, tagValues, dof ) if (ierr .ne. 0) then write(logunit,*) subname,' error in setting ifrac on ice moab instance ' call shr_sys_abort(subname//' ERROR in setting ifrac on ice moab instance ') @@ -931,7 +1047,7 @@ subroutine seq_frac_set(infodata, ice, ocn, fractions_a, fractions_i, fractions_ tagname = 'ofrac'//C_NULL_CHAR tagValues = fractions_i%rAttr(ko,:) - ierr = iMOAB_SetDoubleTagStorageWithGid ( mbixid, tagname, lSize , ent_type, tagValues, GlobalIds ) + ierr = iMOAB_SetDoubleTagStorageWithGid ( mbixid, tagname, lSize , ent_type, tagValues, dof ) if (ierr .ne. 0) then write(logunit,*) subname,' error in setting ofrac on ice moab instance ' call shr_sys_abort(subname//' ERROR in setting ofrac on ice moab instance ') diff --git a/driver-moab/main/seq_map_mod.F90 b/driver-moab/main/seq_map_mod.F90 index 55a974cd0f09..5be18429f7f0 100644 --- a/driver-moab/main/seq_map_mod.F90 +++ b/driver-moab/main/seq_map_mod.F90 @@ -235,7 +235,7 @@ subroutine moab_map_init_rcfile( mbsrc, mbtgt, mbintx, discretization_type, & call shr_sys_abort(subname//' ERROR in loading map file') endif if (seq_comm_iamroot(CPLID)) then - write(logunit,'(2A,I6,4A)') subname,'Result: iMOAB map app ID, maptype, mapfile = ', & + write(logunit,'(2A,I12,4A)') subname,'Result: iMOAB map app ID, maptype, mapfile = ', & mbintx,' ',trim(maptype),' ',trim(mapfile), ', identifier: ', trim(sol_identifier) call shr_sys_flush(logunit) endif @@ -328,6 +328,7 @@ subroutine seq_map_map( mapper, av_s, av_d, fldlist, norm, avwts_s, avwtsfld_s, use iMOAB, only: iMOAB_GetMeshInfo, iMOAB_GetDoubleTagStorage, iMOAB_SetDoubleTagStorage, & iMOAB_GetIntTagStorage, iMOAB_ApplyScalarProjectionWeights, & iMOAB_SendElementTag, iMOAB_ReceiveElementTag, iMOAB_FreeSenderBuffers + use seq_comm_mct, only : num_moab_exports implicit none !----------------------------------------------------- @@ -438,7 +439,7 @@ subroutine seq_map_map( mapper, av_s, av_d, fldlist, norm, avwts_s, avwtsfld_s, #ifdef MOABDEBUG if (seq_comm_iamroot(CPLID)) then write(logunit,*) subname, 'iMOAB mapper ',trim(mapper%mbname), ' iMOAB_mapper nfields', & - nfields, ' fldlist_moab=', trim(fldlist_moab) + nfields, ' fldlist_moab=', trim(fldlist_moab), ' moab step ', num_moab_exports call shr_sys_flush(logunit) endif #endif @@ -496,7 +497,7 @@ subroutine seq_map_map( mapper, av_s, av_d, fldlist, norm, avwts_s, avwtsfld_s, #ifdef MOABDEBUG if (seq_comm_iamroot(CPLID)) then write(logunit, *) subname,' iMOAB mapper rearrange or copy ', mapper%mbname, ' send/recv tags ', trim(fldlist_moab), & - ' mbpresent=', mbpresent, ' mbnorm=', mbnorm + ' mbpresent=', mbpresent, ' mbnorm=', mbnorm, ' moab step:', num_moab_exports call shr_sys_flush(logunit) endif #endif @@ -550,7 +551,8 @@ subroutine seq_map_map( mapper, av_s, av_d, fldlist, norm, avwts_s, avwtsfld_s, endif #ifdef MOABDEBUG if (seq_comm_iamroot(CPLID)) then - write(logunit, *) subname,' iMOAB mapper ', mapper%mbname, ' set norm8wt 1 on source with app id: ', mapper%src_mbid + write(logunit, *) subname,' iMOAB mapper ', mapper%mbname, ' set norm8wt 1 on source with app id: ', & + mapper%src_mbid, ' moab step:', num_moab_exports call shr_sys_flush(logunit) endif #endif @@ -584,7 +586,8 @@ subroutine seq_map_map( mapper, av_s, av_d, fldlist, norm, avwts_s, avwtsfld_s, #ifdef MOABDEBUG if (seq_comm_iamroot(CPLID)) then write(logunit, *) subname,' iMOAB projection mapper: ', mapper%mbname, ' normalize nfields=', & - nfields, ' arrsize_src on root:', arrsize_src, ' shape(targtags_ini)=', shape(targtags_ini) + nfields, ' arrsize_src on root:', arrsize_src, ' shape(targtags_ini)=', shape(targtags_ini), & + ' moab step:', num_moab_exports call shr_sys_flush(logunit) endif #endif @@ -613,7 +616,7 @@ subroutine seq_map_map( mapper, av_s, av_d, fldlist, norm, avwts_s, avwtsfld_s, #ifdef MOABDEBUG if (seq_comm_iamroot(CPLID)) then write(logunit, *) subname,' iMOAB mapper receiving tags with intx and intx_mbid: ', & - mapper%mbname, trim(fldlist_moab) + mapper%mbname, trim(fldlist_moab), ' moab step:', num_moab_exports endif #endif ierr = iMOAB_ReceiveElementTag( mapper%intx_mbid, fldlist_moab, mapper%mpicom, mapper%src_context ) @@ -634,7 +637,8 @@ subroutine seq_map_map( mapper, av_s, av_d, fldlist, norm, avwts_s, avwtsfld_s, #ifdef MOABDEBUG if (seq_comm_iamroot(CPLID)) then - write(logunit, *) subname,' iMOAB projection mapper: ',trim(mapper%mbname), ' between ', mapper%src_mbid, ' and ', mapper%tgt_mbid, trim(fldlist_moab) + write(logunit, *) subname,' iMOAB projection mapper: ',trim(mapper%mbname), ' between ', mapper%src_mbid, ' and ', mapper%tgt_mbid, trim(fldlist_moab), & + ' moab step:', num_moab_exports call shr_sys_flush(logunit) endif #endif diff --git a/driver-moab/main/seq_rest_mod.F90 b/driver-moab/main/seq_rest_mod.F90 index dafd4d056686..5454918d7d01 100644 --- a/driver-moab/main/seq_rest_mod.F90 +++ b/driver-moab/main/seq_rest_mod.F90 @@ -362,7 +362,7 @@ subroutine seq_rest_read(rest_file, infodata, & end subroutine seq_rest_read -subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al) +subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al, samegrid_lr) use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid ! coupler side instances use iMOAB, only: iMOAB_GetGlobalInfo @@ -373,6 +373,7 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al) character(*) , intent(in) :: rest_file ! restart file path/name type(seq_infodata_type), intent(in) :: infodata logical , intent(in) :: samegrid_al ! needed for land nx + logical , intent(in) :: samegrid_lr ! needed for land nx, too integer(IN) :: n,n1,n2,n3 real(r8),allocatable :: ds(:) ! for reshaping diag data for restart file @@ -454,7 +455,12 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al) ierr = iMOAB_GetGlobalInfo(mbaxid, dummy, nx_lnd) ! max id for land will come from atm call seq_io_read(moab_rest_file, mblxid, 'fractions_lx', & 'afrac:lfrac:lfrin', nx=nx_lnd) - else + else if(samegrid_lr) then + ! nx for land will be from global nb rof + ierr = iMOAB_GetGlobalInfo(mbrxid, dummy, nx_lnd) ! max id for land will come from rof + call seq_io_read(moab_rest_file, mblxid, 'fractions_lx', & + 'afrac:lfrac:lfrin', nx=nx_lnd) + else ! is this ever true ? call seq_io_read(moab_rest_file, mblxid, 'fractions_lx', & 'afrac:lfrac:lfrin') endif @@ -471,7 +477,13 @@ subroutine seq_rest_mb_read(rest_file, infodata, samegrid_al) call seq_io_read(moab_rest_file, mblxid, 'l2racc_lx', & trim(tagname), & matrix = p_l2racc_lm, nx=nx_lnd) - else + else if(samegrid_lr) then + ! nx for land will be from global nb rof + ierr = iMOAB_GetGlobalInfo(mbrxid, dummy, nx_lnd) ! max id for land will come from rof + call seq_io_read(moab_rest_file, mblxid, 'l2racc_lx', & + trim(tagname), & + matrix = p_l2racc_lm, nx=nx_lnd) + else call seq_io_read(moab_rest_file, mblxid, 'l2racc_lx', & trim(tagname), & matrix = p_l2racc_lm ) @@ -942,7 +954,7 @@ end subroutine seq_rest_write subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & atm, lnd, ice, ocn, rof, glc, wav, esp, iac, & - tag, samegrid_al, rest_file) + tag, samegrid_al, samegrid_lr, rest_file) use seq_comm_mct, only: mbaxid, mbixid, mboxid, mblxid, mbrxid, mbofxid ! coupler side instances use iMOAB, only: iMOAB_GetGlobalInfo @@ -965,6 +977,7 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & character(len=*) , intent(in) :: tag logical , intent(in) :: samegrid_al ! needed for land nx + logical , intent(in) :: samegrid_lr ! needed for land nx too, for trigrid case character(len=CL) , intent(out) :: rest_file ! Restart filename @@ -1175,6 +1188,12 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & call seq_io_write(rest_file, mblxid, 'fractions_lx', & 'afrac:lfrac:lfrin', & ! seq_frac_mod: character(*),parameter :: fraclist_l = 'afrac:lfrac:lfrin' whead=whead, wdata=wdata, nx=nx_lnd) + else if(samegrid_lr) then + ! nx for land will be from global nb atmosphere + ierr = iMOAB_GetGlobalInfo(mbrxid, dummy, nx_lnd) ! max id for land will come from rof + call seq_io_write(rest_file, mblxid, 'fractions_lx', & + 'afrac:lfrac:lfrin', & ! seq_frac_mod: character(*),parameter :: fraclist_l = 'afrac:lfrac:lfrin' + whead=whead, wdata=wdata, nx=nx_lnd) else call seq_io_write(rest_file, mblxid, 'fractions_lx', & 'afrac:lfrac:lfrin', & ! seq_frac_mod: character(*),parameter :: fraclist_l = 'afrac:lfrac:lfrin' @@ -1197,6 +1216,12 @@ subroutine seq_rest_mb_write(EClock_d, seq_SyncClock, infodata, & call seq_io_write(rest_file, mblxid, 'l2racc_lx', & trim(tagname), & whead=whead, wdata=wdata, matrix = p_l2racc_lm, nx=nx_lnd) + else if(samegrid_lr) then + ! nx for land will be from global nb atmosphere + ierr = iMOAB_GetGlobalInfo(mbrxid, dummy, nx_lnd) ! max id for land will come from rof + call seq_io_write(rest_file, mblxid, 'l2racc_lx', & + trim(tagname), & + whead=whead, wdata=wdata, matrix = p_l2racc_lm, nx=nx_lnd) else call seq_io_write(rest_file, mblxid, 'l2racc_lx', & trim(tagname), & diff --git a/driver-moab/shr/CMakeLists.txt b/driver-moab/shr/CMakeLists.txt index 08d47cd358ce..ebd39b5f4a6e 100644 --- a/driver-moab/shr/CMakeLists.txt +++ b/driver-moab/shr/CMakeLists.txt @@ -1,5 +1,6 @@ list(APPEND drv_sources glc_elevclass_mod.F90 + glc_zocnclass_mod.F90 seq_cdata_mod.F90 seq_comm_mct.F90 seq_infodata_mod.F90 diff --git a/driver-moab/shr/glc_zocnclass_mod.F90 b/driver-moab/shr/glc_zocnclass_mod.F90 new file mode 100644 index 000000000000..50bb4d100528 --- /dev/null +++ b/driver-moab/shr/glc_zocnclass_mod.F90 @@ -0,0 +1,343 @@ +module glc_zocnclass_mod + + !--------------------------------------------------------------------- + ! + ! Purpose: + ! + ! This module contains data and routines for operating on GLC ocean z-level classes. + +#include "shr_assert.h" + use shr_kind_mod, only : r8 => shr_kind_r8 + use shr_sys_mod + use seq_comm_mct, only : logunit + use shr_log_mod, only : errMsg => shr_log_errMsg + + implicit none + save + private + + !-------------------------------------------------------------------------- + ! Public interfaces + !-------------------------------------------------------------------------- + + public :: glc_zocnclass_init ! initialize GLC z-ocean class data + public :: glc_zocnclass_clean ! deallocate memory allocated here + public :: glc_get_num_zocn_classes ! get the number of z-ocean classes + public :: glc_get_zlevels ! get an array of the z-ocean levels + public :: glc_get_zocnclass_bounds ! get the boundaries of all z-ocean classes + public :: glc_zocnclass_as_string ! returns a string corresponding to a given z-ocean class + public :: glc_all_zocnclass_strings ! returns an array of strings for all z-ocean classes + public :: glc_zocn_errcode_to_string ! convert an error code into a string describing the error + + interface glc_zocnclass_init + module procedure glc_zocnclass_init_default + module procedure glc_zocnclass_init_override + end interface glc_zocnclass_init + + + !-------------------------------------------------------------------------- + ! Public data + !-------------------------------------------------------------------------- + + ! Possible error code values + integer, parameter, public :: GLC_ZOCNCLASS_ERR_NONE = 0 ! err_code indicating no error + integer, parameter, public :: GLC_ZOCNCLASS_ERR_UNDEFINED = 1 ! err_code indicating z-ocean classes have not been defined + integer, parameter, public :: GLC_ZOCNCLASS_ERR_TOO_LOW = 2 ! err_code indicating z-level below lowest z-ocean class + integer, parameter, public :: GLC_ZOCNCLASS_ERR_TOO_HIGH = 3 ! err_code indicating z-level above highest z-ocean class + + ! String length for glc z-ocean classes represented as strings + integer, parameter, public :: GLC_ZOCNCLASS_STRLEN = 2 + + !-------------------------------------------------------------------------- + ! Private data + !-------------------------------------------------------------------------- + + ! number of elevation classes + integer :: glc_nzoc ! number of z-ocean classes + + ! z-level of each class. Units are meters above sea level, so values should be <0 + ! indexing goes from shallowest to deepest levels + real(r8), allocatable :: zocn_levels(:) + ! upper and lower z-level limit for each class (m) + ! first dimension: indexing goes from shallowest to deepest levels + ! second dimension: index 1 is upper limit, index 2 is lower limit + real(r8), allocatable :: zocn_bnds(:,:) + + +contains + + !----------------------------------------------------------------------- + subroutine glc_zocnclass_init_default(my_glc_nzoc) + ! + ! !DESCRIPTION: + ! Initialize GLC z-ocean class data to default values, based on given glc_nzoc + ! + ! !USES: + ! + ! !ARGUMENTS: + integer, intent(in) :: my_glc_nzoc ! number of GLC z-ocean classes + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'glc_zocnclass_init' + integer :: i + !----------------------------------------------------------------------- + + glc_nzoc = my_glc_nzoc + allocate(zocn_levels(glc_nzoc)) + + select case (glc_nzoc) + case(0) + ! do nothing + case(4) + zocn_levels = [-250._r8, -750._r8, -1250._r8, -1750._r8] + case(30) + zocn_levels = [ -30._r8, -90._r8, -150._r8, -210._r8, -270._r8, & + -330._r8, -390._r8, -450._r8, -510._r8, -570._r8, & + -630._r8, -690._r8, -750._r8, -810._r8, -870._r8, & + -930._r8, -990._r8, -1050._r8, -1110._r8, -1170._r8, & + -1230._r8, -1290._r8, -1350._r8, -1410._r8, -1470._r8, & + -1530._r8, -1590._r8, -1650._r8, -1710._r8, -1770._r8] + case default + write(logunit,*) subname,' ERROR: unknown glc_nzoc: ', glc_nzoc + call shr_sys_abort(subname//' ERROR: unknown glc_nzoc') + end select + + call glc_zocnclass_init_bnds() + + end subroutine glc_zocnclass_init_default + + !----------------------------------------------------------------------- + subroutine glc_zocnclass_init_bnds() + integer :: i + + allocate(zocn_bnds(2,glc_nzoc)) + zocn_bnds(:,:) = 0._r8 + if (glc_nzoc >= 2) then + zocn_bnds(1,1) = 0._r8 + zocn_bnds(2,1) = 0.5_r8 * (zocn_levels(1) + zocn_levels(2)) + do i = 2, glc_nzoc - 1 + zocn_bnds(1,i) = 0.5_r8 * (zocn_levels(i-1) + zocn_levels(i)) + zocn_bnds(2,i) = 0.5_r8 * (zocn_levels(i) + zocn_levels(i+1)) + enddo + zocn_bnds(1,glc_nzoc) = 0.5_r8 * (zocn_levels(glc_nzoc-1) + zocn_levels(glc_nzoc)) + zocn_bnds(2,glc_nzoc) = zocn_levels(glc_nzoc) + (zocn_levels(glc_nzoc) - zocn_bnds(1,glc_nzoc)) + endif + end subroutine glc_zocnclass_init_bnds + + !----------------------------------------------------------------------- + subroutine glc_zocnclass_init_override(my_glc_nzoc, my_zocn_levels) + ! + ! !DESCRIPTION: + ! Initialize GLC zocn class data to the given z-values + ! + ! The input, my_zocn_levels, should have my_glc_nzoc elements. + ! + ! !USES: + ! + ! !ARGUMENTS: + integer, intent(in) :: my_glc_nzoc ! number of GLC z-ocean classes + real(r8), intent(in) :: my_zocn_levels(:) ! z-ocean values (m) + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'glc_zocnlass_init_override' + !----------------------------------------------------------------------- + + SHR_ASSERT_ALL_FL((ubound(my_zocn_levels) == (/my_glc_nzoc/)), __FILE__, __LINE__) + + glc_nzoc = my_glc_nzoc + allocate(zocn_levels(glc_nzoc)) + zocn_levels = my_zocn_levels + allocate(zocn_bnds(2,glc_nzoc)) + + call glc_zocnclass_init_bnds() + + end subroutine glc_zocnclass_init_override + + !----------------------------------------------------------------------- + subroutine glc_zocnclass_clean() + ! + ! !DESCRIPTION: + ! Deallocate memory allocated in this module + + character(len=*), parameter :: subname = 'glc_zocnclass_clean' + !----------------------------------------------------------------------- + + if (allocated(zocn_levels)) then + deallocate(zocn_levels) + end if + if (allocated(zocn_bnds)) then + deallocate(zocn_bnds) + end if + glc_nzoc = 0 + + end subroutine glc_zocnclass_clean + + !----------------------------------------------------------------------- + function glc_get_num_zocn_classes() result(num_zocn_classes) + ! + ! !DESCRIPTION: + ! Get the number of GLC z-ocean classes + ! + ! !ARGUMENTS: + integer :: num_zocn_classes ! function result + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'glc_get_num_zocn_classes' + !----------------------------------------------------------------------- + + num_zocn_classes = glc_nzoc + + end function glc_get_num_zocn_classes + + !----------------------------------------------------------------------- + function glc_get_zlevels() result(zlevs) + ! + ! !DESCRIPTION: + ! Get all z-levels + ! + ! This returns an array of size (glc_nzoc) + ! + ! !USES: + ! + ! !ARGUMENTS: + real(r8) :: zlevs(glc_nzoc) ! function result + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'glc_get_zlevels' + !----------------------------------------------------------------------- + + zlevs(:) = zocn_levels(:) + + end function glc_get_zlevels + + !----------------------------------------------------------------------- + function glc_get_zocnclass_bounds() result(zocnclass_bounds) + ! + ! !DESCRIPTION: + ! Get the boundaries of all z-ocean classes. + ! + ! This returns an array of size (glc_nzoc,2) + ! + ! !USES: + ! + ! !ARGUMENTS: + real(r8) :: zocnclass_bounds(2,glc_nzoc) ! function result + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'glc_get_zocnclass_bounds' + !----------------------------------------------------------------------- + + zocnclass_bounds(:,:) = zocn_bnds(:,:) + + end function glc_get_zocnclass_bounds + + !----------------------------------------------------------------------- + function glc_zocnclass_as_string(zocn_class) result(zc_string) + ! + ! !DESCRIPTION: + ! Returns a string corresponding to a given elevation class. + ! + ! This string can be used as a suffix for fields in MCT attribute vectors. + ! + ! !USES: + ! + ! !ARGUMENTS: + character(len=GLC_ZOCNCLASS_STRLEN) :: zc_string ! function result + integer, intent(in) :: zocn_class + ! + ! !LOCAL VARIABLES: + character(len=16) :: format_string + + character(len=*), parameter :: subname = 'glc_zocnclass_as_string' + !----------------------------------------------------------------------- + + ! e.g., for GLC_ZOCNCLASS_STRLEN = 2, format_string will be '(i2.2)' + write(format_string,'(a,i0,a,i0,a)') '(i', GLC_ZOCNCLASS_STRLEN, '.', GLC_ZOCNCLASS_STRLEN, ')' + + write(zc_string,trim(format_string)) zocn_class + end function glc_zocnclass_as_string + + !----------------------------------------------------------------------- + function glc_all_zocnclass_strings(include_zero) result(zc_strings) + ! + ! !DESCRIPTION: + ! Returns an array of strings corresponding to all z-ocean classes from 1 to glc_nzoc + ! + ! If include_zero is present and true, then includes z-ocean class 0 - so goes from + ! 0 to glc_nzoc + ! + ! These strings can be used as suffixes for fields in MCT attribute vectors. + ! + ! !USES: + ! + ! !ARGUMENTS: + character(len=GLC_ZOCNCLASS_STRLEN), allocatable :: zc_strings(:) ! function result + logical, intent(in), optional :: include_zero ! if present and true, include elevation class 0 (default is false) + ! + ! !LOCAL VARIABLES: + logical :: l_include_zero ! local version of optional include_zero argument + integer :: lower_bound + integer :: i + + character(len=*), parameter :: subname = 'glc_all_zocnclass_strings' + !----------------------------------------------------------------------- + + if (present(include_zero)) then + l_include_zero = include_zero + else + l_include_zero = .false. + end if + + if (l_include_zero) then + lower_bound = 0 + else + lower_bound = 1 + end if + + allocate(zc_strings(lower_bound:glc_nzoc)) + do i = lower_bound, glc_nzoc + zc_strings(i) = glc_zocnclass_as_string(i) + end do + + end function glc_all_zocnclass_strings + + + !----------------------------------------------------------------------- + function glc_zocn_errcode_to_string(err_code) result(err_string) + ! + ! !DESCRIPTION: + ! + ! + ! !USES: + ! + ! !ARGUMENTS: + character(len=256) :: err_string ! function result + integer, intent(in) :: err_code ! error code (one of the GLC_ZOCNCLASS_ERR* values) + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'glc_errcode_to_string' + !----------------------------------------------------------------------- + + select case (err_code) + case (GLC_ZOCNCLASS_ERR_NONE) + err_string = '(no error)' + case (GLC_ZOCNCLASS_ERR_UNDEFINED) + err_string = 'Z-ocean classes have not yet been defined' + case (GLC_ZOCNCLASS_ERR_TOO_LOW) + err_string = 'Z-level below the lower bound of the lowest z-ocean class' + case (GLC_ZOCNCLASS_ERR_TOO_HIGH) + err_string = 'Z-level above the upper bound of the highest z-ocean class' + case default + err_string = 'UNKNOWN ERROR' + end select + + end function glc_zocn_errcode_to_string + + +end module glc_zocnclass_mod diff --git a/driver-moab/shr/seq_comm_mct.F90 b/driver-moab/shr/seq_comm_mct.F90 index 4d151e61ebca..d434ff3937f8 100644 --- a/driver-moab/shr/seq_comm_mct.F90 +++ b/driver-moab/shr/seq_comm_mct.F90 @@ -238,9 +238,6 @@ module seq_comm_mct integer, public :: mbintxia ! iMOAB id for intx mesh between ice and atmosphere integer, public :: mrofid ! iMOAB id of moab rof app integer, public :: mbrxid ! iMOAB id of moab rof read from file on coupler pes -! integer, public :: mbrmapro ! iMOAB id for read map between river and ocean; it exists on coupler PEs -! ! similar to intx id, oa, la; -! integer, public :: mbrxoid ! iMOAB id for rof migrated to coupler for ocean context (r2o mapping) integer, public :: mbintxro ! iMOAB id for read map between river and ocean; it exists on coupler PEs logical, public :: mbrof_data = .false. ! made true if no rtm mesh, which means data rof ? integer, public :: mbintxar ! iMOAB id for intx mesh between atm and river @@ -678,8 +675,6 @@ subroutine seq_comm_init(global_comm_in, driver_comm_in, nmlfile, drv_comm_id) mbintxia = -1 ! iMOAB id for ice intx with atm on coupler pes mrofid = -1 ! iMOAB id of moab rof app mbrxid = -1 ! iMOAB id of moab rof migrated to coupler - ! mbrmapro = -1 ! iMOAB id of moab instance of map read from rof2ocn map file - ! mbrxoid = -1 ! iMOAB id of moab instance rof to coupler in ocean context mbintxro = -1 ! iMOAB id of moab instance of map read from rof2ocn map file mbintxar = -1 ! iMOAB id for intx mesh between atm and river mbintxlr = -1 ! iMOAB id for intx mesh between land and river diff --git a/driver-moab/shr/seq_flds_mod.F90 b/driver-moab/shr/seq_flds_mod.F90 index bc432558851d..73126149597f 100644 --- a/driver-moab/shr/seq_flds_mod.F90 +++ b/driver-moab/shr/seq_flds_mod.F90 @@ -309,6 +309,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) use shr_string_mod, only : shr_string_listIntersect use shr_mpi_mod, only : shr_mpi_bcast use glc_elevclass_mod, only : glc_elevclass_init + use glc_zocnclass_mod, only : glc_zocnclass_init use seq_infodata_mod, only : seq_infodata_type, seq_infodata_getdata ! !INPUT/OUTPUT PARAMETERS: @@ -404,10 +405,11 @@ subroutine seq_flds_set(nmlfile, ID, infodata) logical :: flds_polar logical :: flds_tf integer :: glc_nec + integer :: glc_nzoc namelist /seq_cplflds_inparm/ & flds_co2a, flds_co2b, flds_co2c, flds_co2_dmsa, flds_wiso, flds_polar, flds_tf, & - glc_nec, ice_ncat, seq_flds_i2o_per_cat, flds_bgc_oi, & + glc_nec, glc_nzoc, ice_ncat, seq_flds_i2o_per_cat, flds_bgc_oi, & nan_check_component_fields, rof_heat, atm_flux_method, atm_gustiness, & rof2ocn_nutrients, lnd_rof_two_way, ocn_rof_two_way, rof_sed, & wav_ocn_coup @@ -447,6 +449,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) flds_polar = .false. flds_tf = .false. glc_nec = 0 + glc_nzoc = 0 ice_ncat = 1 seq_flds_i2o_per_cat = .false. nan_check_component_fields = .false. @@ -483,6 +486,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) call shr_mpi_bcast(flds_polar , mpicom) call shr_mpi_bcast(flds_tf , mpicom) call shr_mpi_bcast(glc_nec , mpicom) + call shr_mpi_bcast(glc_nzoc , mpicom) call shr_mpi_bcast(ice_ncat , mpicom) call shr_mpi_bcast(seq_flds_i2o_per_cat, mpicom) call shr_mpi_bcast(nan_check_component_fields, mpicom) @@ -496,6 +500,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) call shr_mpi_bcast(wav_ocn_coup, mpicom) call glc_elevclass_init(glc_nec) + call glc_zocnclass_init(glc_nzoc) !--------------------------------------------------------------------------- ! Read in namelists for user specified new fields @@ -3014,18 +3019,33 @@ subroutine seq_flds_set(nmlfile, ID, infodata) attname = 'So_rhoeff' call metadata_set(attname, longname, stdname, units) - if (flds_tf) then + if ((flds_tf) .and. (glc_nzoc > 0)) then + ! glc fields with multiple ocn z classes: ocn->glc + ! + ! Note that these fields are sent in multiple elevation classes from ocn->cpl + ! and cpl->ocn (which differs from glc_nec variables) - name = 'So_tf2d' - call seq_flds_add(o2x_states,trim(name)) - call seq_flds_add(x2g_states,trim(name)) - call seq_flds_add(x2g_tf_states_from_ocn,trim(name)) - longname = 'ocean thermal forcing at predefined critical depth' - stdname = 'ocean_thermal_forcing_at_critical_depth' + name = 'So_tf3d' + longname = 'ocean thermal forcing at z-level' + stdname = 'ocean_thermal_forcing_at_z_level' units = 'C' attname = name - call metadata_set(attname, longname, stdname, units) - + call set_glc_zocnclass_field(name, attname, longname, stdname, units, o2x_states) + call set_glc_zocnclass_field(name, attname, longname, stdname, units, x2g_states, & + additional_list = .true.) + call set_glc_zocnclass_field(name, attname, longname, stdname, units, x2g_tf_states_from_ocn, & + additional_list = .true.) + + name = 'So_tf3d_mask' + longname = 'mask of valid ocean thermal forcing at z-level' + stdname = 'mask_ocean_thermal_forcing_at_z_level' + units = 'none' + attname = name + call set_glc_zocnclass_field(name, attname, longname, stdname, units, o2x_states) + call set_glc_zocnclass_field(name, attname, longname, stdname, units, x2g_states, & + additional_list = .true.) + call set_glc_zocnclass_field(name, attname, longname, stdname, units, x2g_tf_states_from_ocn, & + additional_list = .true.) end if name = 'Fogx_qicelo' @@ -4393,6 +4413,66 @@ end subroutine set_glc_elevclass_field !=============================================================================== + subroutine set_glc_zocnclass_field(name, attname, longname, stdname, units, fieldlist, & + additional_list) + + ! Sets a coupling field for all ocn z classes (1:glc_nzoc) + ! + ! Note that, if glc_nzoc = 0, then we don't create any coupling fields + ! + ! Puts the coupling fields in the given fieldlist, and also does the appropriate + ! metadata_set calls. + ! + ! additional_list should be .false. (or absent) the first time this is called for a + ! given set of coupling fields. However, if this same set of coupling fields is being + ! added to multiple field lists, then additional_list should be set to true for the + ! second and subsequent calls; in this case, the metadata_set calls are not done + ! (because they have already been done). + ! + ! name, attname and longname give the base name of the field; the ocn z class + ! index will be appended as a suffix + + ! !USES: + use glc_zocnclass_mod, only : glc_get_num_zocn_classes, glc_zocnclass_as_string + + ! !INPUT/OUTPUT PARAMETERS: + character(len=*), intent(in) :: name ! base field name to add to fieldlist + character(len=*), intent(in) :: attname ! base field name for metadata + character(len=*), intent(in) :: longname ! base long name for metadata + character(len=*), intent(in) :: stdname ! standard name for metadata + character(len=*), intent(in) :: units ! units for metadata + character(len=*), intent(inout) :: fieldlist ! field list into which the fields should be added + + logical, intent(in), optional :: additional_list ! whether this is an additional list for the same set of coupling fields (see above for details; defaults to false) + + !EOP + integer :: num + character(len= 16) :: cnum + logical :: l_additional_list ! local version of the optional additional_list argument + + l_additional_list = .false. + if (present(additional_list)) then + l_additional_list = additional_list + end if + + if (glc_get_num_zocn_classes() > 0) then + do num = 1, glc_get_num_zocn_classes() + cnum = glc_zocnclass_as_string(num) + + call seq_flds_add(fieldlist, trim(name) // trim(cnum)) + + if (.not. l_additional_list) then + call metadata_set(attname = trim(attname) // trim(cnum), & + longname = trim(longname) // ' of thermal forcing class ' // trim(cnum), & + stdname = stdname, & + units = units) + end if + end do + end if + end subroutine set_glc_zocnclass_field + + !=============================================================================== + subroutine seq_flds_esmf_metadata_get(shortname, longname, stdname, units) ! !USES: diff --git a/driver-moab/shr/seq_infodata_mod.F90 b/driver-moab/shr/seq_infodata_mod.F90 index a602ad9b53d3..f8d1f2cb909a 100644 --- a/driver-moab/shr/seq_infodata_mod.F90 +++ b/driver-moab/shr/seq_infodata_mod.F90 @@ -209,6 +209,7 @@ MODULE seq_infodata_mod logical :: glcice_present ! does glc have iceberg coupling on logical :: glc_prognostic ! does component model need input data from driver logical :: glc_coupled_fluxes ! does glc send fluxes to other components (only relevant if glc_present is .true.) + integer :: glc_nzoc ! number of z-levels for ocn/glc thermal forcing coupling logical :: wav_present ! does component model exist logical :: wav_prognostic ! does component model need input data from driver logical :: esp_present ! does component model exist @@ -777,6 +778,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid, cpl_tag) ! if glc_present is .false., so it's okay to just start out assuming it's .true. ! in all cases. infodata%glc_coupled_fluxes = .true. + infodata%glc_nzoc = 0 infodata%wav_prognostic = .false. infodata%iac_prognostic = .false. infodata%iceberg_prognostic = .false. @@ -1024,7 +1026,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ ice_present, ice_prognostic, & glc_present, glc_prognostic, & iac_present, iac_prognostic, & - glc_coupled_fluxes, & + glc_coupled_fluxes, glc_nzoc, & flood_present, wav_present, wav_prognostic, rofice_present, & glclnd_present, glcocn_present, glcice_present, iceberg_prognostic,& esp_present, esp_prognostic, & @@ -1205,6 +1207,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ logical, optional, intent(OUT) :: glcice_present logical, optional, intent(OUT) :: glc_prognostic logical, optional, intent(OUT) :: glc_coupled_fluxes + integer, optional, intent(OUT) :: glc_nzoc logical, optional, intent(OUT) :: wav_present logical, optional, intent(OUT) :: wav_prognostic logical, optional, intent(OUT) :: iac_present @@ -1399,6 +1402,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ if ( present(glcice_present) ) glcice_present = infodata%glcice_present if ( present(glc_prognostic) ) glc_prognostic = infodata%glc_prognostic if ( present(glc_coupled_fluxes)) glc_coupled_fluxes = infodata%glc_coupled_fluxes + if ( present(glc_nzoc) ) glc_nzoc = infodata%glc_nzoc if ( present(wav_present) ) wav_present = infodata%wav_present if ( present(wav_prognostic) ) wav_prognostic = infodata%wav_prognostic if ( present(esp_present) ) esp_present = infodata%esp_present @@ -1593,7 +1597,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ ocn_c2_glcshelf, ocn_c2_glctf, & ice_present, ice_prognostic, & glc_present, glc_prognostic, & - glc_coupled_fluxes, & + glc_coupled_fluxes, glc_nzoc, & flood_present, wav_present, wav_prognostic, rofice_present, & glclnd_present, glcocn_present, glcice_present, iceberg_prognostic,& esp_present, esp_prognostic, & @@ -1776,6 +1780,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ logical, optional, intent(IN) :: glcice_present logical, optional, intent(IN) :: glc_prognostic logical, optional, intent(IN) :: glc_coupled_fluxes + integer, optional, intent(IN) :: glc_nzoc logical, optional, intent(IN) :: wav_present logical, optional, intent(IN) :: wav_prognostic logical, optional, intent(IN) :: esp_present @@ -1969,6 +1974,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ if ( present(glcice_present) ) infodata%glcice_present = glcice_present if ( present(glc_prognostic) ) infodata%glc_prognostic = glc_prognostic if ( present(glc_coupled_fluxes)) infodata%glc_coupled_fluxes = glc_coupled_fluxes + if ( present(glc_nzoc) ) infodata%glc_nzoc = glc_nzoc if ( present(wav_present) ) infodata%wav_present = wav_present if ( present(wav_prognostic) ) infodata%wav_prognostic = wav_prognostic if ( present(iac_present) ) infodata%iac_present = iac_present @@ -2285,6 +2291,7 @@ subroutine seq_infodata_bcast(infodata,mpicom) call shr_mpi_bcast(infodata%glcice_present, mpicom) call shr_mpi_bcast(infodata%glc_prognostic, mpicom) call shr_mpi_bcast(infodata%glc_coupled_fluxes, mpicom) + call shr_mpi_bcast(infodata%glc_nzoc, mpicom) call shr_mpi_bcast(infodata%wav_present, mpicom) call shr_mpi_bcast(infodata%wav_prognostic, mpicom) call shr_mpi_bcast(infodata%esp_present, mpicom) @@ -2602,6 +2609,7 @@ subroutine seq_infodata_Exchange(infodata,ID,type) call shr_mpi_bcast(infodata%glcice_present, mpicom, pebcast=cmppe) call shr_mpi_bcast(infodata%glc_prognostic, mpicom, pebcast=cmppe) call shr_mpi_bcast(infodata%glc_coupled_fluxes, mpicom, pebcast=cmppe) + call shr_mpi_bcast(infodata%glc_nzoc, mpicom, pebcast=cmppe) call shr_mpi_bcast(infodata%glc_nx, mpicom, pebcast=cmppe) call shr_mpi_bcast(infodata%glc_ny, mpicom, pebcast=cmppe) ! dead_comps is true if it's ever set to true @@ -2661,6 +2669,7 @@ subroutine seq_infodata_Exchange(infodata,ID,type) call shr_mpi_bcast(infodata%glcice_present, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%glc_prognostic, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%glc_coupled_fluxes, mpicom, pebcast=cplpe) + call shr_mpi_bcast(infodata%glc_nzoc, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%wav_present, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%wav_prognostic, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%iac_present, mpicom, pebcast=cplpe) @@ -3015,6 +3024,7 @@ SUBROUTINE seq_infodata_print( infodata ) write(logunit,F0L) subname,'glcice_present = ', infodata%glcice_present write(logunit,F0L) subname,'glc_prognostic = ', infodata%glc_prognostic write(logunit,F0L) subname,'glc_coupled_fluxes = ', infodata%glc_coupled_fluxes + write(logunit,F0L) subname,'glc_nzoc = ', infodata%glc_nzoc write(logunit,F0L) subname,'wav_present = ', infodata%wav_present write(logunit,F0L) subname,'wav_prognostic = ', infodata%wav_prognostic write(logunit,F0L) subname,'iac_present = ', infodata%iac_present