diff --git a/components/elm/src/cpl/lnd_comp_mct.F90 b/components/elm/src/cpl/lnd_comp_mct.F90 index e533d5bcbc2b..8d3ae5e2997b 100644 --- a/components/elm/src/cpl/lnd_comp_mct.F90 +++ b/components/elm/src/cpl/lnd_comp_mct.F90 @@ -18,7 +18,6 @@ module lnd_comp_mct #ifdef HAVE_MOAB use seq_comm_mct, only: mlnid! id of moab land app - use seq_comm_mct, only: mb_land_mesh! true if land is full mesh (on the river mesh) use seq_comm_mct, only: num_moab_exports #ifdef MOABCOMP use seq_comm_mct , only: seq_comm_compare_mb_mct @@ -51,7 +50,6 @@ module lnd_comp_mct integer :: mpicom_lnd_moab ! used also for mpi-reducing the difference between moab tags and mct avs integer :: rank2 - logical :: samegrid_al ! #endif !--------------------------------------------------------------------------- @@ -314,13 +312,6 @@ subroutine lnd_init_mct( EClock, cdata_l, x2l_l, l2x_l, NLFilename ) call lnd_domain_mct( bounds, lsz, gsMap_lnd, dom_l ) #ifdef HAVE_MOAB -! find out samegrid_al or not; from infodata - samegrid_al = .true. - call seq_infodata_GetData(infodata , & - atm_gnam=atm_gnam , & - lnd_gnam=lnd_gnam ) - if (trim(atm_gnam) /= trim(lnd_gnam)) samegrid_al = .false. - mb_land_mesh = .not. samegrid_al ! global variable, saved in seq_comm call init_moab_land(bounds, LNDID) #endif call mct_aVect_init(x2l_l, rList=seq_flds_x2l_fields, lsize=lsz) @@ -547,8 +538,7 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l) ! loop over all fields in seq_flds_x2l_fields call mct_list_init(temp_list ,seq_flds_x2l_fields) size_list=mct_list_nitem (temp_list) - ent_type = 0 ! entity type is vertex for land, usually (bigrid case) - if (mb_land_mesh) ent_type = 1 + ent_type = 0 ! entity type is vertex for land, always if (rank2 .eq. 0) print *, num_moab_exports, trim(seq_flds_x2l_fields), ' lnd import check' do index_list = 1, size_list call mct_list_get(mctOStr,index_list,temp_list) @@ -846,7 +836,7 @@ subroutine init_moab_land(bounds, LNDID) use spmdmod , only: masterproc use iMOAB , only: iMOAB_CreateVertices, iMOAB_WriteMesh, iMOAB_RegisterApplication, & iMOAB_DefineTagStorage, iMOAB_SetIntTagStorage, iMOAB_SetDoubleTagStorage, & - iMOAB_ResolveSharedEntities, iMOAB_CreateElements, iMOAB_UpdateMeshInfo + iMOAB_ResolveSharedEntities, iMOAB_UpdateMeshInfo type(bounds_type) , intent(in) :: bounds integer , intent(in) :: LNDID ! id of the land app @@ -893,200 +883,97 @@ subroutine init_moab_land(bounds, LNDID) vgids(n) = ldecomp%gdc2glo(bounds%begg+n-1) ! local to global ! end do gsize = ldomain%ni * ldomain%nj ! size of the total grid - ! if ldomain%nv > 3 , create mesh - - ! Case where land and river share mesh (tri-grid) - if (ldomain%nv .ge. 3 .and. .not.samegrid_al) then - ! number of vertices is nv * lsz ! - allocate(moab_vert_coords(lsz*dims*ldomain%nv)) - ! loop over ldomain - allocate(moabconn(ldomain%nv * lsz)) - do n = bounds%begg, bounds%endg - i = (n - bounds%begg) * ldomain%nv - do iv = 1, ldomain%nv - lonv = ldomain%mblonv(n, iv) * SHR_CONST_PI/180. - latv = ldomain%mblatv(n, iv) * SHR_CONST_PI/180. - - i = i + 1 ! iv-th vertex of cell n; i starts at 1 - moab_vert_coords(3*i-2)=COS(latv)*COS(lonv) - moab_vert_coords(3*i-1)=COS(latv)*SIN(lonv) - moab_vert_coords(3*i )=SIN(latv) - moabconn(i) = i - enddo - enddo - ierr = iMOAB_CreateVertices(mlnid, lsz * 3 * ldomain%nv, dims, moab_vert_coords) - if (ierr > 0 ) & - call endrun('Error: fail to create MOAB vertices in land model') - - mbtype = 2 ! triangle - if (ldomain%nv .eq. 4) mbtype = 3 ! quad - if (ldomain%nv .gt. 4) mbtype = 4 ! polygon - block_ID = 100 !some value - ierr = iMOAB_CreateElements( mlnid, lsz, mbtype, ldomain%nv, moabconn, block_ID ); - - - ! define some useful tags on cells - tagtype = 0 ! dense, integer - numco = 1 - tagname='GLOBAL_ID'//C_NULL_CHAR - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to retrieve GLOBAL_ID tag ') - - ent_type = 1 ! element type - ierr = iMOAB_SetIntTagStorage ( mlnid, tagname, lsz , ent_type, vgids) - if (ierr > 0 ) & - call endrun('Error: fail to set GLOBAL_ID tag ') - - ! use moab_vert_coords as a data holder for a frac tag and area tag that we will create - ! on the vertices; do not allocate other data array - ! Define and Set Fraction - tagname='frac'//C_NULL_CHAR - tagtype = 1 ! dense, double - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to create frac tag ') - - do i = 1, lsz - n = i-1 + bounds%begg - moab_vert_coords(i) = ldomain%frac(n) - enddo - ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz , ent_type, moab_vert_coords) - if (ierr > 0 ) & - call endrun('Error: fail to set frac tag ') - - ! Define and Set area - tagname='area'//C_NULL_CHAR - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to create area tag ') - do i = 1, lsz - n = i-1 + bounds%begg - moab_vert_coords(i) = ldomain%area(n)/(re*re) ! use the same doubles for second tag :) - enddo - - ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz , ent_type, moab_vert_coords ) - if (ierr > 0 ) & - call endrun('Error: fail to set area tag ') - - ! Define aream - tagname='aream'//C_NULL_CHAR - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to create aream tag ') - - deallocate(moabconn) - deallocate(vgids) - - - ! Now do the verticies - allocate(vgids(lsz*ldomain%nv)) ! - do n = 1, lsz - do i=1,ldomain%nv - vgids( (n-1)*ldomain%nv+i ) = (ldecomp%gdc2glo(bounds%begg+n-1)-1)*ldomain%nv+i ! local to global ! - end do - end do - ent_type = 0 ! vertices now - tagname = 'GLOBAL_ID'//C_NULL_CHAR - ierr = iMOAB_SetIntTagStorage ( mlnid, tagname, lsz , ent_type, vgids ) - if (ierr > 0 ) & - call endrun('Error: fail to set global ID tag on vertices in land mesh ') - ierr = iMOAB_UpdateMeshInfo( mlnid ) - if (ierr > 0 ) & - call endrun('Error: fail to update mesh info ') - - ! Case where land and atmosphere share mesh - else ! old point cloud mesh - allocate(moab_vert_coords(lsz*dims)) - do i = 1, lsz - n = i-1 + bounds%begg - lonv = ldomain%lonc(n) *SHR_CONST_PI/180. - latv = ldomain%latc(n) *SHR_CONST_PI/180. - moab_vert_coords(3*i-2)=COS(latv)*COS(lonv) - moab_vert_coords(3*i-1)=COS(latv)*SIN(lonv) - moab_vert_coords(3*i )=SIN(latv) - enddo - ierr = iMOAB_CreateVertices(mlnid, lsz*3, dims, moab_vert_coords) - if (ierr > 0 ) & - call endrun('Error: fail to create MOAB vertices in land model') - - tagtype = 0 ! dense, integer - numco = 1 - tagname='GLOBAL_ID'//C_NULL_CHAR - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to retrieve GLOBAL_ID tag ') - - ent_type = 0 ! vertex type - ierr = iMOAB_SetIntTagStorage ( mlnid, tagname, lsz , ent_type, vgids) - if (ierr > 0 ) & - call endrun('Error: fail to set GLOBAL_ID tag ') - - ierr = iMOAB_ResolveSharedEntities( mlnid, lsz, vgids ); - if (ierr > 0 ) & - call endrun('Error: fail to resolve shared entities') - - !there are no shared entities, but we will set a special partition tag, in order to see the - ! partitions ; it will be visible with a Pseudocolor plot in VisIt - tagname='partition'//C_NULL_CHAR - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to create new partition tag ') - - vgids = iam - ierr = iMOAB_SetIntTagStorage ( mlnid, tagname, lsz , ent_type, vgids) - if (ierr > 0 ) & - call endrun('Error: fail to set partition tag ') - - ! use moab_vert_coords as a data holder for a frac tag and area tag that we will create - ! on the vertices; do not allocate other data array - tagname='frac'//C_NULL_CHAR - tagtype = 1 ! dense, double - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to create frac tag ') - - do i = 1, lsz - n = i-1 + bounds%begg - moab_vert_coords(i) = ldomain%frac(n) - enddo - ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz , ent_type, moab_vert_coords) - if (ierr > 0 ) & - call endrun('Error: fail to set frac tag ') - - tagname='area'//C_NULL_CHAR - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to create area tag ') - do i = 1, lsz - n = i-1 + bounds%begg - moab_vert_coords(i) = ldomain%area(n)/(re*re) ! use the same doubles for second tag :) - enddo - - ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz , ent_type, moab_vert_coords ) - if (ierr > 0 ) & - call endrun('Error: fail to set area tag ') - - ! aream needed in cime_init for now. - tagname='aream'//C_NULL_CHAR - ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) - if (ierr > 0 ) & - call endrun('Error: fail to create aream tag ') - ! ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz , ent_type, moab_vert_coords ) - ! if (ierr > 0 ) & - ! call endrun('Error: fail to set aream tag ') - ierr = iMOAB_UpdateMeshInfo( mlnid ) - if (ierr > 0 ) & - call endrun('Error: fail to update mesh info ') - endif - ! add more domain fields that are missing from domain fields: lat, lon, mask, hgt - tagname = 'lat:lon:mask:hgt'//C_NULL_CHAR - tagtype = 1 ! dense, double + + allocate(moab_vert_coords(lsz*dims)) + do i = 1, lsz + n = i-1 + bounds%begg + lonv = ldomain%lonc(n) *SHR_CONST_PI/180. + latv = ldomain%latc(n) *SHR_CONST_PI/180. + moab_vert_coords(3*i-2)=COS(latv)*COS(lonv) + moab_vert_coords(3*i-1)=COS(latv)*SIN(lonv) + moab_vert_coords(3*i )=SIN(latv) + enddo + ierr = iMOAB_CreateVertices(mlnid, lsz*3, dims, moab_vert_coords) + if (ierr > 0 ) & + call endrun('Error: fail to create MOAB vertices in land model') + + tagtype = 0 ! dense, integer numco = 1 + tagname='GLOBAL_ID'//C_NULL_CHAR + ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call endrun('Error: fail to retrieve GLOBAL_ID tag ') + + ent_type = 0 ! vertex type + ierr = iMOAB_SetIntTagStorage ( mlnid, tagname, lsz , ent_type, vgids) + if (ierr > 0 ) & + call endrun('Error: fail to set GLOBAL_ID tag ') + + ierr = iMOAB_ResolveSharedEntities( mlnid, lsz, vgids ); + if (ierr > 0 ) & + call endrun('Error: fail to resolve shared entities') + + !there are no shared entities, but we will set a special partition tag, in order to see the + ! partitions ; it will be visible with a Pseudocolor plot in VisIt + tagname='partition'//C_NULL_CHAR + ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call endrun('Error: fail to create new partition tag ') + + vgids = iam + ierr = iMOAB_SetIntTagStorage ( mlnid, tagname, lsz , ent_type, vgids) + if (ierr > 0 ) & + call endrun('Error: fail to set partition tag ') + + ! use moab_vert_coords as a data holder for a frac tag and area tag that we will create + ! on the vertices; do not allocate other data array + tagname='frac'//C_NULL_CHAR + tagtype = 1 ! dense, double + ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call endrun('Error: fail to create frac tag ') + + do i = 1, lsz + n = i-1 + bounds%begg + moab_vert_coords(i) = ldomain%frac(n) + enddo + ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz , ent_type, moab_vert_coords) + if (ierr > 0 ) & + call endrun('Error: fail to set frac tag ') + + tagname='area'//C_NULL_CHAR ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) if (ierr > 0 ) & - call endrun('Error: fail to create lat:lon:mask:hgt tags ') + call endrun('Error: fail to create area tag ') + do i = 1, lsz + n = i-1 + bounds%begg + moab_vert_coords(i) = ldomain%area(n)/(re*re) ! use the same doubles for second tag :) + enddo + + ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz , ent_type, moab_vert_coords ) + if (ierr > 0 ) & + call endrun('Error: fail to set area tag ') + + ! aream needed in cime_init for now. + tagname='aream'//C_NULL_CHAR + ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call endrun('Error: fail to create aream tag ') + ! ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz , ent_type, moab_vert_coords ) + ! if (ierr > 0 ) & + ! call endrun('Error: fail to set aream tag ') + ierr = iMOAB_UpdateMeshInfo( mlnid ) + if (ierr > 0 ) & + call endrun('Error: fail to update mesh info ') + + ! add more domain fields that are missing from domain fields: lat, lon, mask, hgt + tagname = 'lat:lon:mask:hgt'//C_NULL_CHAR + tagtype = 1 ! dense, double + numco = 1 + ierr = iMOAB_DefineTagStorage(mlnid, tagname, tagtype, numco, tagindex ) + if (ierr > 0 ) & + call endrun('Error: fail to create lat:lon:mask:hgt tags ') ! moab_vert_coords is big enough in both case to hold enough data for us: lat, lon, mask do i = 1, lsz @@ -1098,9 +985,7 @@ subroutine init_moab_land(bounds, LNDID) tagname = 'lat:lon:mask'//C_NULL_CHAR ent_type = 0 ! point cloud usually - if (ldomain%nv .ge. 3 .and. .not.samegrid_al) then - ent_type = 1 ! cell in tri-grid case - endif + ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsz*3 , ent_type, moab_vert_coords) if (ierr > 0 ) & call endrun('Error: fail to set lat lon mask tag ') @@ -1255,11 +1140,7 @@ subroutine lnd_export_moab(EClock, bounds, lnd2atm_vars, lnd2glc_vars) end do tagname=trim(seq_flds_l2x_fields)//C_NULL_CHAR - if (samegrid_al) then - ent_type = 0 ! vertices, cells only if samegrid_al false - else - ent_type = 1 - endif + ent_type = 0 ! vertices only, from now on ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, totalmbls , ent_type, l2x_lm(1,1) ) if (ierr > 0 ) & call shr_sys_abort( sub//' Error: fail to set moab l2x '// trim(seq_flds_l2x_fields) ) @@ -1452,11 +1333,7 @@ subroutine lnd_import_moab(EClock, bounds, atm2lnd_vars, glc2lnd_vars) call endrun('Error: fail to write the moab lnd mesh before import ') #endif tagname=trim(seq_flds_x2l_fields)//C_NULL_CHAR - if (samegrid_al) then - ent_type = 0 ! vertices, cells only if samegrid_al false - else - ent_type = 1 - endif + ent_type = 0 ! vertices ierr = iMOAB_GetDoubleTagStorage ( mlnid, tagname, totalmblsimp , ent_type, x2l_lm(1,1) ) if ( ierr > 0) then call endrun('Error: fail to get seq_flds_x2l_fields for land moab instance on component') diff --git a/components/elm/src/main/surfrdMod.F90 b/components/elm/src/main/surfrdMod.F90 index a8146e5a0f15..97a4fae62d48 100755 --- a/components/elm/src/main/surfrdMod.F90 +++ b/components/elm/src/main/surfrdMod.F90 @@ -20,11 +20,6 @@ module surfrdMod use ncdio_pio , only : ncd_io, check_var, ncd_inqfdims, check_dim, ncd_inqdid, ncd_inqdlen use pio -#ifdef HAVE_MOAB - use mct_mod , only : mct_gsMap - use decompMod , only : get_elmlevel_gsmap - ! use spmdMod , only : iam ! rank on the land communicator -#endif use spmdMod use topounit_varcon , only : max_topounits, has_topounit @@ -184,11 +179,6 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename) ! pflotran:beg----------------------------- integer :: j, np, nv -#ifdef HAVE_MOAB - type(mct_gsMap), pointer :: gsMap - integer :: i, iv , iseg, ig, local ! ni, nj, nv, nseg, global ig - -#endif ! pflotran:end----------------------------- character(len=32) :: subname = 'surfrd_get_grid' ! subroutine name @@ -258,59 +248,6 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename) end if ! pflotran:end----------------------------------------------- - -#ifdef HAVE_MOAB - ! read xv and yv for MOAB to learn mesh verticies - if (ldomain%nv>=3 ) then - call get_elmlevel_gsmap (grlnd, gsMap) - allocate(rdata3d(nv,ni,nj)) ! transpose from c, as this is fortran - vname = 'xv' - ! this should be improved in a distributed read, that does not use full grid ni * nj * nv 720*360*4*8 ~ 8Mb - call ncd_io(ncid=ncid, varname=trim(vname), data=rdata3d, flag='read', readvar=readvar) - if (.not. readvar) call endrun( msg=trim(subname)//' ERROR: xv NOT on file'//errMsg(__FILE__, __LINE__)) - ! fill up the ldomain%mblonv(begg:endg, 1:nv) array - local = begg - do iseg = 1, gsMap%ngseg - if (gsMap%pe_loc(iseg) .eq. iam) then - do ig = gsMap%start(iseg), gsMap%start(iseg) + gsMap%length(iseg) - 1 - j = (ig-1)/ni + 1 - i = ig - ni*(j-1) - do iv = 1, nv - if (local .le. endg) then - ldomain%mblonv(local, iv ) = rdata3d(iv, i, j) - else - write (iulog, *), 'OVERFLOW', iseg, gsMap%pe_loc(iseg), gsMap%start(iseg), gsMap%length(iseg), local - endif - enddo - local = local + 1 - enddo - endif - enddo - ! repeat for mblatv - vname = 'yv' - call ncd_io(ncid=ncid, varname=trim(vname), data=rdata3d, flag='read', readvar=readvar) - if (.not. readvar) call endrun( msg=trim(subname)//' ERROR: yv NOT on file'//errMsg(__FILE__, __LINE__)) - ! fill up the ldomain%lonv(begg:endg, 1:nv) array - local = begg - do iseg = 1, gsMap%ngseg - if (gsMap%pe_loc(iseg) .eq. iam) then - do ig = gsMap%start(iseg), gsMap%start(iseg) + gsMap%length(iseg) - 1 - j = (ig-1)/ni + 1 - i = ig - ni*(j-1) - do iv = 1, nv - if (local .le. endg) then - ldomain%mblatv(local, iv ) = rdata3d(iv, i, j) - endif - enddo - local = local + 1 - enddo - endif - enddo - ! deallocate what is not needed anymore (for half degree land model, ~8Mb) - deallocate(rdata3d) - - end if -#endif else call ncd_io(ncid=ncid, varname= 'AREA', flag='read', data=ldomain%area, & dim1name=grlnd, readvar=readvar) diff --git a/components/elm/src/utils/domainMod.F90 b/components/elm/src/utils/domainMod.F90 index 5ef3ae611cf6..2c7771179d21 100755 --- a/components/elm/src/utils/domainMod.F90 +++ b/components/elm/src/utils/domainMod.F90 @@ -52,10 +52,6 @@ module domainMod integer :: nv ! number of vertices real(r8),pointer :: latv(:,:) ! latitude of grid cell's vertices (deg) real(r8),pointer :: lonv(:,:) ! longitude of grid cell's vertices (deg) -#ifdef HAVE_MOAB - real(r8),pointer :: mblatv(:,:) ! latitude of grid cell's vertices (deg) for MOAB - real(r8),pointer :: mblonv(:,:) ! longitude of grid cell's vertices (deg) for MOAB -#endif real(r8) :: lon0 ! the origin lon/lat (Most western/southern corner, if not globally covered grids; OR -180W(360E)/-90N) real(r8) :: lat0 ! the origin lon/lat (Most western/southern corner, if not globally covered grids; OR -180W(360E)/-90N) @@ -154,22 +150,6 @@ subroutine domain_init(domain,isgrid2d,ni,nj,nbeg,nend,elmlevel) endif end if ! pflotran:end----------------------------------------------------- -#ifdef HAVE_MOAB - if (domain%nv > 0 .and. domain%nv /= huge(1)) then - if(.not.associated(domain%mblonv)) then - allocate(domain%mblonv(nb:ne, 1:domain%nv), stat=ier) - if (ier /= 0) & - call shr_sys_abort('domain_init ERROR: allocate mblonv ') - domain%mblonv = nan - endif - if(.not.associated(domain%mblatv)) then - allocate(domain%mblatv(nb:ne, 1:domain%nv)) - if (ier /= 0) & - call shr_sys_abort('domain_init ERROR: allocate mblatv ') - domain%mblatv = nan - endif - end if -#endif if (present(elmlevel)) then domain%elmlevel = elmlevel @@ -265,23 +245,6 @@ subroutine domain_clean(domain) endif endif ! pflotran:beg----------------------------------------------------- -#ifdef HAVE_MOAB - if (domain%nv > 0 .and. domain%nv /= huge(1)) then - if (associated(domain%mblonv)) then - deallocate(domain%mblonv, stat=ier) - if (ier /= 0) & - call shr_sys_abort('domain_clean ERROR: deallocate mblonv ') - nullify(domain%mblonv) - endif - - if (associated(domain%mblatv)) then - deallocate(domain%mblatv, stat=ier) - if (ier /= 0) & - call shr_sys_abort('domain_clean ERROR: deallocate mblatv ') - nullify(domain%mblatv) - endif - endif -#endif else if (masterproc) then diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F index d1b140563bb7..f13e2c2e5c70 100644 --- a/components/mpas-ocean/driver/ocn_comp_mct.F +++ b/components/mpas-ocean/driver/ocn_comp_mct.F @@ -884,13 +884,7 @@ end subroutine xml_stream_get_attributes ! !----------------------------------------------------------------------- - call ocn_import_mct(x2o_o, errorCode) - if (errorCode /= 0) then - call mpas_log_write('Error in ocn_import_mct', MPAS_LOG_CRIT) - endif - #ifdef HAVE_MOAB - #ifdef MOABCOMP ! loop over all fields in seq_flds_x2o_fields call mct_list_init(temp_list ,seq_flds_x2o_fields) @@ -906,6 +900,14 @@ end subroutine xml_stream_get_attributes enddo call mct_list_clean(temp_list) #endif +#endif + + call ocn_import_mct(x2o_o, errorCode) + if (errorCode /= 0) then + call mpas_log_write('Error in ocn_import_mct', MPAS_LOG_CRIT) + endif + +#ifdef HAVE_MOAB call ocn_import_moab(Eclock, errorCode) if (errorCode /= 0) then @@ -1029,12 +1031,8 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{ call mpas_get_timeInterval(timeStep, dt=dt) call mpas_reset_clock_alarm(domain_ptr % clock, coupleAlarmID, ierr=ierr) - ! Import state from coupler - call ocn_import_mct(x2o_o, ierr) - ! Import state from moab coupler -#ifdef HAVE_MOAB - +#ifdef HAVE_MOAB #ifdef MOABCOMP ! loop over all fields in seq_flds_x2o_fields call mct_list_init(temp_list ,seq_flds_x2o_fields) @@ -1049,6 +1047,12 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{ enddo call mct_list_clean(temp_list) #endif +#endif + + ! Import state from coupler + call ocn_import_mct(x2o_o, ierr) + ! Import state from moab coupler +#ifdef HAVE_MOAB call ocn_import_moab(Eclock, ierr) if (ierr /= 0) then @@ -3146,15 +3150,19 @@ end subroutine datetime!}}} #ifdef HAVE_MOAB -! import method from moab -! copied from ocn_import_mct, will replace x2o_o AV with x2o_om array read locally - subroutine ocn_import_moab( Eclock, errorCode)!{{{ + +!*********************************************************************** +!BOP +! !IROUTINE: ocn_import_moab +! !INTERFACE: + + subroutine ocn_import_moab(Eclock, errorCode)!{{{ ! !DESCRIPTION: !----------------------------------------------------------------------- -! This routine receives message from cpl7 driver +! This routine receives message from moab driver ! -! The following fields are always received from the coupler: +! The following fields are always received from the driver: ! ! o taux -- zonal wind stress (taux) (W/m2 ) ! o tauy -- meridonal wind stress (tauy) (W/m2 ) @@ -3171,11 +3179,22 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ ! o ifrac -- ice fraction (%) ! o rofl -- river runoff flux (kg/m2/s) ! o rofi -- ice runoff flux (kg/m2/s) +! o rofDIN -- DIN runoff flux (kg/m2/s) +! o rofDIP -- DIP runoff flux (kg/m2/s) +! o rofDON -- DON runoff flux (kg/m2/s) +! o rofDOP -- DOP runoff flux (kg/m2/s) +! o rofDOC -- DOC runoff flux (kg/m2/s) +! o rofPP -- PP runoff flux (kg/m2/s) +! o rofDSi -- DSi runoff flux (kg/m2/s) +! o rofPOC -- POC runoff flux (kg/m2/s) +! o rofPN -- PN runoff flux (kg/m2/s) +! o rofDIC -- DIC runoff flux (kg/m2/s) +! o rofFe -- Fe runoff flux (kg/m2/s) ! ! The following fields are sometimes received from the coupler, ! depending on model options: ! -! o pbot -- bottom atm pressure (Pa) +! o pslv -- atmospheric pressure at sea level (Pa) ! o duu10n -- 10m wind speed squared (m^2/s^2) ! o co2prog-- bottom atm level prognostic co2 ! o co2diag-- bottom atm level diagnostic co2 @@ -3184,28 +3203,15 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ ! ! !REVISION HISTORY: ! same as module - + use iMOAB, only : iMOAB_GetDoubleTagStorage, iMOAB_WriteMesh ! !INPUT/OUTPUT PARAMETERS: + type(ESMF_Clock), intent(inout) :: EClock ! type(mct_aVect) , intent(inout) :: x2o_o ! instead, we will get x2o_om from MPOID ! !OUTPUT PARAMETERS: - use iMOAB, only : iMOAB_GetDoubleTagStorage, iMOAB_WriteMesh - !EOP - !BOC - !----------------------------------------------------------------------- - ! - ! local variables - !----------------------------------------------------------------------- - ! - ! local variables - ! - !----------------------------------------------------------------------- - integer :: ent_type, ierr - character(CXX) :: tagname - type(ESMF_Clock), intent(inout) :: EClock integer, intent(out) :: & errorCode ! returned error code @@ -3221,6 +3227,9 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ label, & message + integer :: ent_type, ierr + character(CXX) :: tagname + integer :: & i,n @@ -3232,6 +3241,7 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ config_use_DMSTracers_sea_ice_coupling, & config_use_MacroMoleculesTracers, & config_use_MacroMoleculesTracers_sea_ice_coupling, & + config_use_CFCTracers, & config_remove_ais_river_runoff, & config_remove_ais_ice_runoff, & config_cvmix_kpp_use_theory_wave @@ -3250,7 +3260,8 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ ecosysAuxiliary, & ecosysSeaIceCoupling, & DMSSeaIceCoupling, & - MacroMoleculesSeaIceCoupling + MacroMoleculesSeaIceCoupling, & + CFCAuxiliary integer, pointer :: nCellsSolve @@ -3278,10 +3289,22 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ iceFluxFeParticulateField, & iceFluxFeDissolvedField, & iceFluxDustField, & + riverFluxNO3Field, & + riverFluxPO4Field, & + riverFluxSiO3Field, & + riverFluxDOCField, & + riverFluxDONField, & + riverFluxDOPField, & + riverFluxDICField, & + riverFluxALKField, & + riverFluxFeField, & landIceFreshwaterFluxField, & landIceHeatFluxField, & landIceFractionField, & - windSpeed10mField + windSpeed10mField, & + significantWaveHeightField, & + peakWaveFrequencyField, & + peakWaveDirectionField !landIcePressureField type (field2DReal), pointer :: iceFluxPhytoCField, & @@ -3289,6 +3312,9 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ type (field2DReal), pointer :: landIceInterfaceTracersField + type (field2DReal), pointer :: stokesDriftZonalWavenumberField, & + stokesDriftMeridionalWavenumberField + real (kind=RKIND), dimension(:), pointer :: windStressZonal, windStressMeridional, & latentHeatFlux, sensibleHeatFlux, & longWaveHeatFluxUp, & @@ -3302,6 +3328,7 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ atmosphericPressure, iceFraction, & seaIcePressure, windSpeedSquared10m, & atmosphericCO2, atmosphericCO2_ALT_CO2, & + windSpeedSquared10mCFC, & iceFluxDIC, & iceFluxDON, & iceFluxNO3, & @@ -3313,26 +3340,44 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ iceFluxFeParticulate, & iceFluxFeDissolved, & iceFluxDust, & + riverFluxNO3, & + riverFluxPO4, & + riverFluxSiO3, & + riverFluxDOC, & + riverFluxDON, & + riverFluxDOP, & + riverFluxDIC, & + riverFluxALK, & + riverFluxFe, & landIceFreshwaterFlux, & landIceHeatFlux, & landIceFraction, & - windSpeed10m + areaCell, & + windSpeed10m, & + significantWaveHeight, & + peakWaveFrequency, & + peakWaveDirection !landIcePressure real (kind=RKIND), dimension(:), pointer :: latCell real (kind=RKIND), dimension(:,:), pointer :: iceFluxPhytoC, & - iceFluxDOC + iceFluxDOC, & + stokesDriftZonalWavenumber, & + stokesDriftMeridionalWavenumber real (kind=RKIND) :: removedRiverRunoffFluxThisProc, removedIceRunoffFluxThisProc real (kind=RKIND) :: removedRiverRunoffFluxReduced, removedIceRunoffFluxReduced real (kind=RKIND), dimension(:,:), pointer :: landIceInterfaceTracers + real (kind=RKIND) :: riverFactor + !----------------------------------------------------------------------- ! ! zero out padded cells ! +!----------------------------------------------------------------------- !----------------------------------------------------------------------- integer :: cur_ocn_stepno #ifdef MOABDEBUG @@ -3349,11 +3394,10 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ write(ocnLogUnit,*) 'Fail to write ocean state ' endif #endif + errorCode = 0 ! get moab tags from MPOID - - ent_type = 1 ! cells ! get all tags in one method tagname = trim(seq_flds_x2o_fields)//C_NULL_CHAR @@ -3361,7 +3405,6 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ if ( ierr /= 0 ) then write(ocnLogUnit,*) 'Fail to get MOAB fields ' endif - !----------------------------------------------------------------------- ! ! unpack and distribute wind stress, then convert to correct units @@ -3380,6 +3423,7 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ config_use_DMSTracers_sea_ice_coupling) call mpas_pool_get_config(domain % configs, 'config_use_MacroMoleculesTracers_sea_ice_coupling', & config_use_MacroMoleculesTracers_sea_ice_coupling) + call mpas_pool_get_config(domain % configs, 'config_use_CFCTracers', config_use_CFCTracers) call mpas_pool_get_config(domain % configs, 'config_remove_ais_river_runoff', config_remove_ais_river_runoff) call mpas_pool_get_config(domain % configs, 'config_remove_ais_ice_runoff', config_remove_ais_ice_runoff) call mpas_pool_get_config(domain % configs, 'config_cvmix_kpp_use_theory_wave', config_cvmix_kpp_use_theory_wave) @@ -3418,6 +3462,11 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ call mpas_pool_get_field(forcingPool, 'iceRunoffFlux', iceRunoffFluxField) call mpas_pool_get_field(forcingPool, 'removedRiverRunoffFlux', removedRiverRunoffFluxField) call mpas_pool_get_field(forcingPool, 'removedIceRunoffFlux', removedIceRunoffFluxField) + call mpas_pool_get_field(forcingPool, 'stokesDriftZonalWavenumber', stokesDriftZonalWavenumberField) + call mpas_pool_get_field(forcingPool, 'stokesDriftMeridionalWavenumber', stokesDriftMeridionalWavenumberField) + call mpas_pool_get_field(forcingPool, 'significantWaveHeight', significantWaveHeightField) + call mpas_pool_get_field(forcingPool, 'peakWaveFrequency', peakWaveFrequencyField) + call mpas_pool_get_field(forcingPool, 'peakWaveDirection', peakWaveDirectionField) call mpas_pool_get_field(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFluxField) call mpas_pool_get_field(forcingPool, 'landIceHeatFlux', landIceHeatFluxField) @@ -3459,9 +3508,15 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ landIceInterfaceTracers => landIceInterfaceTracersField % array landIceFraction => landIceFractionField % array windSpeed10m => windSpeed10mField % array + stokesDriftZonalWavenumber => stokesDriftZonalWavenumberField % array + stokesDriftMeridionalWavenumber => stokesDriftMeridionalWavenumberField % array + significantWaveHeight => significantWaveHeightField % array + peakWaveFrequency => peakWaveFrequencyField % array + peakWaveDirection => peakWaveDirectionField % array !landIcePressure => landIcePressureField % array call mpas_pool_get_array(meshPool, 'latCell', latCell) + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) ! BGC fields if (config_use_ecosysTracers) then @@ -3474,6 +3529,27 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ call mpas_pool_get_field(ecosysAuxiliary, 'atmosphericCO2_ALT_CO2', atmosphericCO2_ALT_CO2Field) atmosphericCO2_ALT_CO2 => atmosphericCO2_ALT_CO2Field % array + if (config_use_ecosysTracers_river_inputs_from_coupler) then + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxNO3' , riverFluxNO3Field) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxPO4' , riverFluxPO4Field) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxDON' , riverFluxDONField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxDOP' , riverFluxDOPField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxSiO3', riverFluxSiO3Field) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxDOC' , riverFluxDOCField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxDIC' , riverFluxDICField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxALK' , riverFluxALKField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxFe' , riverFluxFeField) + riverFluxNO3 => riverFluxNO3Field % array + riverFluxPO4 => riverFluxPO4Field % array + riverFluxDON => riverFluxDONField % array + riverFluxDOP => riverFluxDOPField % array + riverFluxSiO3 => riverFluxSiO3Field % array + riverFluxDOC => riverFluxDOCField % array + riverFluxDIC => riverFluxDICField % array + riverFluxALK => riverFluxALKField % array + riverFluxFe => riverFluxFeField % array + endif + call mpas_pool_get_config(domain % configs, 'config_ecosys_atm_co2_option', & config_ecosys_atm_co2_option) call mpas_pool_get_config(domain % configs, 'config_ecosys_atm_alt_co2_option', & @@ -3519,6 +3595,13 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ iceFluxDMSP => iceFluxDMSPField % array endif + ! CFC fields + if (config_use_CFCTracers) then + call mpas_pool_get_subpool(forcingPool, 'CFCAuxiliary', CFCAuxiliary) + call mpas_pool_get_field(CFCAuxiliary, 'windSpeedSquared10mCFC', windSpeedSquared10mField) + windSpeedSquared10mCFC => windSpeedSquared10mField % array + endif + if (config_remove_ais_river_runoff) then ! Initialize this field removedRiverRunoffFlux(:) = 0.0_RKIND @@ -3533,8 +3616,7 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ ! Initialize this field windSpeed10m(:) = 0.0_RKIND endif - -! replace 'x2o_o % rAttr(' to 'x2o_om(n, ' and ', n)' with ')' +! ! replace 'x2o_o % rAttr(' to 'x2o_om(n, ' and ', n)' with ')' do i = 1, nCellsSolve n = n + 1 if ( windStressZonalField % isActive ) then @@ -3618,6 +3700,32 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ iceFraction(i) = x2o_om(n, index_x2o_Si_ifrac) end if + if ( stokesDriftZonalWavenumberField % isActive ) then + stokesDriftZonalWavenumber(1,i) = x2o_om(n, index_x2o_Sw_ustokes_wavenumber_1) + stokesDriftZonalWavenumber(2,i) = x2o_om(n, index_x2o_Sw_ustokes_wavenumber_2) + stokesDriftZonalWavenumber(3,i) = x2o_om(n, index_x2o_Sw_ustokes_wavenumber_3) + stokesDriftZonalWavenumber(4,i) = x2o_om(n, index_x2o_Sw_ustokes_wavenumber_4) + stokesDriftZonalWavenumber(5,i) = x2o_om(n, index_x2o_Sw_ustokes_wavenumber_5) + stokesDriftZonalWavenumber(6,i) = x2o_om(n, index_x2o_Sw_ustokes_wavenumber_6) + end if + if ( stokesDriftMeridionalWavenumberField % isActive ) then + stokesDriftMeridionalWavenumber(1,i) = x2o_om(n, index_x2o_Sw_vstokes_wavenumber_1) + stokesDriftMeridionalWavenumber(2,i) = x2o_om(n, index_x2o_Sw_vstokes_wavenumber_2) + stokesDriftMeridionalWavenumber(3,i) = x2o_om(n, index_x2o_Sw_vstokes_wavenumber_3) + stokesDriftMeridionalWavenumber(4,i) = x2o_om(n, index_x2o_Sw_vstokes_wavenumber_4) + stokesDriftMeridionalWavenumber(5,i) = x2o_om(n, index_x2o_Sw_vstokes_wavenumber_5) + stokesDriftMeridionalWavenumber(6,i) = x2o_om(n, index_x2o_Sw_vstokes_wavenumber_6) + end if + if ( significantWaveHeightField % isActive ) then + significantWaveHeight(i) = x2o_om(n, index_x2o_Sw_Hs) + end if + if ( peakWaveFrequencyField % isActive ) then + peakWaveFrequency(i) = x2o_om(n, index_x2o_Sw_Fp) + end if + if ( peakWaveDirectionField % isActive ) then + peakWaveDirection(i) = x2o_om(n, index_x2o_Sw_Dp) + end if + if (config_cvmix_kpp_use_theory_wave) then if ( windSpeed10mField% isActive ) then windSpeed10m(i) = sqrt( x2o_om(n, index_x2o_So_duu10n)) @@ -3675,12 +3783,36 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ else if ( config_ecosys_atm_alt_co2_option == 'bdrc') then atmosphericCO2_ALT_CO2(i) = config_ecosys_atm_co2_constant_value else if ( config_ecosys_atm_alt_co2_option == 'bdrd') then - atmosphericCO2_ALT_CO2(i) = x2o_om(n, index_x2o_Sa_co2diag) + atmosphericCO2_ALT_CO2(i) = x2o_om(n, index_x2o_Sa_co2diag) else atmosphericCO2_ALT_CO2(i) = config_ecosys_atm_co2_constant_value end if end if + if (config_use_ecosysTracers_river_inputs_from_coupler) then + riverFluxNO3(i) = x2o_om(n, index_x2o_Foxx_rofDIN) + riverFluxPO4(i) = x2o_om(n, index_x2o_Foxx_rofDIP) + riverFluxDON(i) = x2o_om(n, index_x2o_Foxx_rofDON) + riverFluxDOP(i) = x2o_om(n, index_x2o_Foxx_rofDOP) + riverFluxSiO3(i) = x2o_om(n, index_x2o_Foxx_rofDSi) + riverFluxDOC(i) = x2o_om(n, index_x2o_Foxx_rofDOC) + riverFluxDIC(i) = x2o_om(n, index_x2o_Foxx_rofDIC) + riverFluxFe(i) = x2o_om(n, index_x2o_Foxx_rofFe ) + +! convert from kgNutrient/(m2-s) to mmol/m3 m/s + riverFactor = 1.e6_RKIND + riverFluxNO3(i) = riverFluxNO3(i)*riverFactor/14.007_RKIND + riverFluxPO4(i) = riverFluxPO4(i)*riverFactor/30.974_RKIND + riverFluxDON(i) = riverFluxDON(i)*riverFactor/14.007_RKIND + riverFluxDOP(i) = riverFluxDOP(i)*riverFactor/30.974_RKIND + riverFluxSiO3(i) = riverFluxSiO3(i)*riverFactor/28.085_RKIND + riverFluxDOC(i) = riverFluxDOC(i)*riverFactor/12.001_RKIND + riverFluxDIC(i) = riverFluxDIC(i)*riverFactor/12.001_RKIND + riverFluxFe(i) = riverFluxFe(i)*riverFactor/55.845_RKIND + + riverFluxALK(i) = riverFluxDIC(i) + endif + if (config_use_ecosysTracers_sea_ice_coupling) then if ( iceFluxPhytoCField % isActive ) then iceFluxPhytoC(1,i) = x2o_om(n, index_x2o_Fioi_algae1) @@ -3714,6 +3846,7 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ if ( iceFluxDOCField % isActive ) then iceFluxDOC(1,i) = x2o_om(n, index_x2o_Fioi_doc1) iceFluxDOC(2,i) = x2o_om(n, index_x2o_Fioi_doc2) + iceFluxDOC(3,i) = x2o_om(n, index_x2o_Fioi_doc3) endif if ( iceFluxDONField % isActive ) then iceFluxDON(i) = x2o_om(n, index_x2o_Fioi_don1) @@ -3730,6 +3863,13 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ endif endif + ! CFC fields + if (config_use_CFCTracers) then + if ( windSpeedSquared10mField % isActive ) then + windSpeedSquared10mCFC(i) = x2o_om(n, index_x2o_So_duu10n) + end if + end if + end do block_ptr => block_ptr % next @@ -3757,6 +3897,11 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ call mpas_pool_get_field(forcingPool, 'atmosphericPressure', atmosphericPressureField) call mpas_pool_get_field(forcingPool, 'seaIcePressure', seaIcePressureField) call mpas_pool_get_field(forcingPool, 'iceFraction', iceFractionField) + call mpas_pool_get_field(forcingPool, 'stokesDriftZonalWavenumber', stokesDriftZonalWavenumberField) + call mpas_pool_get_field(forcingPool, 'stokesDriftMeridionalWavenumber', stokesDriftMeridionalWavenumberField) + call mpas_pool_get_field(forcingPool, 'significantWaveHeight', significantWaveHeightField) + call mpas_pool_get_field(forcingPool, 'peakWaveFrequency', peakWaveFrequencyField) + call mpas_pool_get_field(forcingPool, 'peakWaveDirection', peakWaveDirectionField) call mpas_pool_get_field(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFluxField) call mpas_pool_get_field(forcingPool, 'landIceHeatFlux', landIceHeatFluxField) @@ -3777,6 +3922,18 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ call mpas_pool_get_field(ecosysAuxiliary, 'atmosphericCO2', atmosphericCO2Field) call mpas_pool_get_field(ecosysAuxiliary, 'atmosphericCO2_ALT_CO2', atmosphericCO2_ALT_CO2Field) + if (config_use_ecosysTracers_river_inputs_from_coupler) then + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxNO3' , riverFluxNO3Field) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxPO4' , riverFluxPO4Field) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxDON' , riverFluxDONField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxDOP' , riverFluxDOPField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxSiO3', riverFluxSiO3Field) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxDOC' , riverFluxDOCField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxDIC' , riverFluxDICField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxALK' , riverFluxALKField) + call mpas_pool_get_field(ecosysAuxiliary, 'riverFluxFe' , riverFluxFeField) + endif + if (config_use_ecosysTracers_sea_ice_coupling) then call mpas_pool_get_subpool(forcingPool, 'ecosysSeaIceCoupling', ecosysSeaIceCoupling) @@ -3800,6 +3957,12 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ call mpas_pool_get_field(DMSSeaIceCoupling, 'iceFluxDMSP', iceFluxDMSPField) endif + ! CFC fields + if (config_use_CFCTracers) then + call mpas_pool_get_subpool(forcingPool, 'CFCAuxiliary', CFCAuxiliary) + call mpas_pool_get_field(CFCAuxiliary, 'windSpeedSquared10mCFC', windSpeedSquared10mField) + endif + if ( windStressMeridionalField % isActive ) then call mpas_dmpar_exch_halo_field(windStressMeridionalField) end if @@ -3860,6 +4023,21 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ if ( iceFractionField % isActive ) then call mpas_dmpar_exch_halo_field(iceFractionField) end if + if ( stokesDriftZonalWavenumberField % isActive ) then + call mpas_dmpar_exch_halo_field(stokesDriftZonalWavenumberField) + end if + if ( stokesDriftMeridionalWavenumberField % isActive ) then + call mpas_dmpar_exch_halo_field(stokesDriftMeridionalWavenumberField) + end if + if ( significantWaveHeightField % isActive ) then + call mpas_dmpar_exch_halo_field(significantWaveHeightField) + end if + if ( peakWaveFrequencyField % isActive ) then + call mpas_dmpar_exch_halo_field(peakWaveFrequencyField) + end if + if ( peakWaveDirectionField % isActive ) then + call mpas_dmpar_exch_halo_field(peakWaveDirectionField) + end if if ( landIceFreshwaterFluxField % isActive ) then call mpas_dmpar_exch_halo_field(landIceFreshwaterFluxField) @@ -3894,6 +4072,36 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ call mpas_dmpar_exch_halo_field(atmosphericCO2_ALT_CO2Field) end if + if (config_use_ecosysTracers_river_inputs_from_coupler) then + if ( riverFluxNO3Field % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxNO3Field) + end if + if ( riverFluxPO4Field % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxPO4Field) + end if + if ( riverFluxDONField % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxDONField) + end if + if ( riverFluxDOPField % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxDOPField) + end if + if ( riverFluxSiO3Field % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxSiO3Field) + end if + if ( riverFluxDOCField % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxDOCField) + end if + if ( riverFluxDICField % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxDICField) + end if + if ( riverFluxALKField % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxALKField) + end if + if ( riverFluxFeField % isActive ) then + call mpas_dmpar_exch_halo_field(riverFluxFeField) + end if + endif + if (config_use_ecosysTracers_sea_ice_coupling) then if ( iceFluxPhytoCField % isActive ) then call mpas_dmpar_exch_halo_field(iceFluxPhytoCField) @@ -3939,302 +4147,348 @@ subroutine ocn_import_moab( Eclock, errorCode)!{{{ endif endif + ! CFC fields + if (config_use_CFCTracers .and. .not. config_use_ecosysTracers) then + if ( windSpeedSquared10mField % isActive ) then + call mpas_dmpar_exch_halo_field(windSpeedSquared10mField) + end if + endif + !----------------------------------------------------------------------- !EOC end subroutine ocn_import_moab!}}} +!*********************************************************************** +!BOP +! !IROUTINE: ocn_export_moab +! !INTERFACE: - subroutine ocn_export_moab(EClock) !{{{ - - ! !DESCRIPTION: - ! This routine calls the routines necessary to send mpas ocean fields to MOAB coupler - ! - use iMOAB, only : iMOAB_SetDoubleTagStorage, iMOAB_WriteMesh - !EOP - !BOC - type(ESMF_Clock) , intent(inout) :: EClock ! Input synchronization clock from driver - ! - ! local variables - ! - !----------------------------------------------------------------------- - integer :: ent_type, ierr, cur_ocn_stepno - character(len=100) :: outfile, wopts, localmeshfile, lnum - character(CXX) :: tagname - - integer :: i, n - integer, pointer :: nCellsSolve, index_temperatureSurfaceValue, index_salinitySurfaceValue, & - index_avgZonalSurfaceVelocity, index_avgMeridionalSurfaceVelocity, & - index_avgZonalSSHGradient, index_avgMeridionalSSHGradient + subroutine ocn_export_moab(EClock) !{{{ - type (block_type), pointer :: block_ptr +! !DESCRIPTION: +! This routine calls the routines necessary to send MPASO fields to +! the MOAB driver +! +! !REVISION HISTORY: +! same as module + use iMOAB, only : iMOAB_SetDoubleTagStorage, iMOAB_WriteMesh +! !INPUT/OUTPUT PARAMETERS: - type (mpas_pool_type), pointer :: meshPool, & - forcingPool, & - statePool, & - tracersPool, & - ecosysAuxiliary, & - ecosysSeaIceCoupling, & - DMSSeaIceCoupling, & - MacroMoleculesSeaIceCoupling - - integer, dimension(:), pointer :: landIceMask - - real (kind=RKIND), dimension(:), pointer :: seaIceEnergy, accumulatedFrazilIceMass, frazilSurfacePressure, & - avgTotalFreshWaterTemperatureFlux, & - avgCO2_gas_flux, DMSFlux, surfaceUpwardCO2Flux, & - avgOceanSurfaceDIC, & - avgOceanSurfaceDON, & - avgOceanSurfaceNO3, & - avgOceanSurfaceSiO3, & - avgOceanSurfaceNH4, & - avgOceanSurfaceDMS, & - avgOceanSurfaceDMSP, & - avgOceanSurfaceDOCr, & - avgOceanSurfaceDOCSemiLabile, & - avgOceanSurfaceFeParticulate, & - avgOceanSurfaceFeDissolved, & - ssh, & - avgLandIceFreshwaterFlux, & - avgRemovedRiverRunoffFlux, & - avgRemovedIceRunoffFlux, & - avgLandIceHeatFlux, & - avgRemovedIceRunoffHeatFlux - - real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, & - avgSSHGradient, avgOceanSurfacePhytoC, & - avgOceanSurfaceDOC, layerThickness - - real (kind=RKIND) :: surfaceFreezingTemp - - logical, pointer :: frazilIceActive, & - config_remove_ais_river_runoff, & - config_remove_ais_ice_runoff, & - config_use_ecosysTracers, & - config_use_DMSTracers, & - config_use_MacroMoleculesTracers, & - config_use_ecosysTracers_sea_ice_coupling, & - config_use_DMSTracers_sea_ice_coupling, & - config_use_MacroMoleculesTracers_sea_ice_coupling - - character (len=StrKIND), pointer :: config_land_ice_flux_mode - - logical :: keepFrazil - - - ! get configure options - call mpas_pool_get_package(domain % packages, 'frazilIceActive', frazilIceActive) - call mpas_pool_get_config(domain % configs, 'config_use_ecosysTracers', config_use_ecosysTracers) - call mpas_pool_get_config(domain % configs, 'config_land_ice_flux_mode', config_land_ice_flux_mode) - call mpas_pool_get_config(domain % configs, 'config_remove_ais_river_runoff', config_remove_ais_river_runoff) - call mpas_pool_get_config(domain % configs, 'config_remove_ais_ice_runoff', config_remove_ais_ice_runoff) - call mpas_pool_get_config(domain % configs, 'config_use_DMSTracers', config_use_DMSTracers) - call mpas_pool_get_config(domain % configs, 'config_use_MacroMoleculesTracers', config_use_MacroMoleculesTracers) - call mpas_pool_get_config(domain % configs, 'config_use_ecosysTracers_sea_ice_coupling', & - config_use_ecosysTracers_sea_ice_coupling) - call mpas_pool_get_config(domain % configs, 'config_use_DMSTracers_sea_ice_coupling', & - config_use_DMSTracers_sea_ice_coupling) - call mpas_pool_get_config(domain % configs, 'config_use_MacroMoleculesTracers_sea_ice_coupling', & - config_use_MacroMoleculesTracers_sea_ice_coupling) - - n = 0 - block_ptr => domain % blocklist - do while(associated(block_ptr)) - call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) - call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) - call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool) + type(ESMF_Clock) , intent(inout) :: EClock ! Input synchronization clock from driver - call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) +! !OUTPUT PARAMETERS: - call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) +!EOP +!BOC +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + integer :: ent_type, ierr, cur_ocn_stepno + character(len=100) :: outfile, wopts, localmeshfile, lnum + character(CXX) :: tagname - call mpas_pool_get_dimension(forcingPool, 'index_avgTemperatureSurfaceValue', index_temperatureSurfaceValue) - call mpas_pool_get_dimension(forcingPool, 'index_avgSalinitySurfaceValue', index_salinitySurfaceValue) - call mpas_pool_get_dimension(forcingPool, 'index_avgSurfaceVelocityZonal', index_avgZonalSurfaceVelocity) - call mpas_pool_get_dimension(forcingPool, 'index_avgSurfaceVelocityMeridional', index_avgMeridionalSurfaceVelocity) - call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientZonal', index_avgZonalSSHGradient) - call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientMeridional', index_avgMeridionalSSHGradient) + integer :: i, n + integer, pointer :: nCellsSolve, index_temperatureSurfaceValue, index_salinitySurfaceValue, & + index_avgZonalSurfaceVelocity, index_avgMeridionalSurfaceVelocity, & + index_avgZonalSSHGradient, index_avgMeridionalSSHGradient + type (block_type), pointer :: block_ptr - call mpas_pool_get_array(statePool, 'ssh', ssh, 1) - call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) + type (mpas_pool_type), pointer :: meshPool, & + forcingPool, & + statePool, & + tracersPool, & + ecosysAuxiliary, & + ecosysSeaIceCoupling, & + DMSSeaIceCoupling, & + MacroMoleculesSeaIceCoupling - call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask) - call mpas_pool_get_array(forcingPool, 'avgTracersSurfaceValue', avgTracersSurfaceValue) - call mpas_pool_get_array(forcingPool, 'avgSurfaceVelocity', avgSurfaceVelocity) - call mpas_pool_get_array(forcingPool, 'avgSSHGradient', avgSSHGradient) - call mpas_pool_get_array(forcingPool, 'avgTotalFreshWaterTemperatureFlux', avgTotalFreshWaterTemperatureFlux) - if ( frazilIceActive ) then - call mpas_pool_get_array(forcingPool, 'seaIceEnergy', seaIceEnergy) - call mpas_pool_get_array(forcingPool, 'frazilSurfacePressure', frazilSurfacePressure) - call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMass, 1) - end if + integer, dimension(:), pointer :: landIceMask - if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then - call mpas_pool_get_array(forcingPool, 'avgLandIceFreshwaterFlux', avgLandIceFreshwaterFlux) - call mpas_pool_get_array(forcingPool, 'avgLandIceHeatFlux', avgLandIceHeatFlux) - endif - if (config_remove_ais_river_runoff) then - call mpas_pool_get_array(forcingPool, 'avgRemovedRiverRunoffFlux', avgRemovedRiverRunoffFlux) - endif - if (config_remove_ais_ice_runoff) then - call mpas_pool_get_array(forcingPool, 'avgRemovedIceRunoffFlux', avgRemovedIceRunoffFlux) - call mpas_pool_get_array(forcingPool, 'avgRemovedIceRunoffHeatFlux', avgRemovedIceRunoffHeatFlux) - endif + real (kind=RKIND), dimension(:), pointer :: seaIceEnergy, accumulatedFrazilIceMass, frazilSurfacePressure, & + avgTotalFreshWaterTemperatureFlux, & + avgCO2_gas_flux, DMSFlux, surfaceUpwardCO2Flux, & + avgOceanSurfaceDIC, & + avgOceanSurfaceDON, & + avgOceanSurfaceNO3, & + avgOceanSurfaceSiO3, & + avgOceanSurfaceNH4, & + avgOceanSurfaceDMS, & + avgOceanSurfaceDMSP, & + avgOceanSurfaceDOCr, & + avgOceanSurfaceDOCSemiLabile, & + avgOceanSurfaceFeParticulate, & + avgOceanSurfaceFeDissolved, & + ssh, & + avgLandIceFreshwaterFlux, & + avgRemovedRiverRunoffFlux, & + avgRemovedIceRunoffFlux, & + avgLandIceHeatFlux, & + avgRemovedIceRunoffHeatFlux - ! BGC fields - if (config_use_ecosysTracers) then + real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, & + avgSSHGradient, avgOceanSurfacePhytoC, & + avgOceanSurfaceDOC, layerThickness - call mpas_pool_get_subpool(forcingPool, 'ecosysAuxiliary', ecosysAuxiliary) - call mpas_pool_get_array(ecosysAuxiliary, 'avgCO2_gas_flux', avgCO2_gas_flux) + real (kind=RKIND) :: surfaceFreezingTemp - end if + logical, pointer :: frazilIceActive, & + config_remove_ais_river_runoff, & + config_remove_ais_ice_runoff, & + config_use_ecosysTracers, & + config_use_DMSTracers, & + config_use_MacroMoleculesTracers, & + config_use_ecosysTracers_sea_ice_coupling, & + config_use_DMSTracers_sea_ice_coupling, & + config_use_MacroMoleculesTracers_sea_ice_coupling - if (config_use_ecosysTracers .and. config_use_ecosysTracers_sea_ice_coupling) then - call mpas_pool_get_subpool(forcingPool, 'ecosysSeaIceCoupling', ecosysSeaIceCoupling) + character (len=StrKIND), pointer :: config_land_ice_flux_mode - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfacePhytoC', avgOceanSurfacePhytoC) - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceDIC', avgOceanSurfaceDIC) - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceNO3', avgOceanSurfaceNO3) - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceSiO3', avgOceanSurfaceSiO3) - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceNH4', avgOceanSurfaceNH4) - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceDOCr', avgOceanSurfaceDOCr) - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceDOCSemiLabile', avgOceanSurfaceDOCSemiLabile) - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceFeParticulate', avgOceanSurfaceFeParticulate) - call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceFeDissolved', avgOceanSurfaceFeDissolved) - endif - if (config_use_DMSTracers .and. config_use_DMSTracers_sea_ice_coupling) then - call mpas_pool_get_subpool(forcingPool, 'DMSSeaIceCoupling', DMSSeaIceCoupling) + logical :: keepFrazil - call mpas_pool_get_array(DMSSeaIceCoupling, 'avgOceanSurfaceDMS', avgOceanSurfaceDMS) - call mpas_pool_get_array(DMSSeaIceCoupling, 'avgOceanSurfaceDMSP', avgOceanSurfaceDMSP) - endif - if (config_use_MacroMoleculesTracers .and. config_use_MacroMoleculesTracers_sea_ice_coupling) then - call mpas_pool_get_subpool(forcingPool, 'MacroMoleculesSeaIceCoupling', MacroMoleculesSeaIceCoupling) - call mpas_pool_get_array(MacroMoleculesSeaIceCoupling, 'avgOceanSurfaceDOC', avgOceanSurfaceDOC) - call mpas_pool_get_array(MacroMoleculesSeaIceCoupling, 'avgOceanSurfaceDON', avgOceanSurfaceDON) - endif - ! call mpas_pool_get_array(forcingPool, 'CO2Flux', CO2Flux) - ! call mpas_pool_get_array(forcingPool, 'DMSFlux', DMSFlux) - ! call mpas_pool_get_array(forcingPool, 'surfaceUpwardCO2Flux', surfaceUpwardCO2Flux) + ! get configure options + call mpas_pool_get_package(domain % packages, 'frazilIceActive', frazilIceActive) + call mpas_pool_get_config(domain % configs, 'config_use_ecosysTracers', config_use_ecosysTracers) + call mpas_pool_get_config(domain % configs, 'config_land_ice_flux_mode', config_land_ice_flux_mode) + call mpas_pool_get_config(domain % configs, 'config_remove_ais_river_runoff', config_remove_ais_river_runoff) + call mpas_pool_get_config(domain % configs, 'config_remove_ais_ice_runoff', config_remove_ais_ice_runoff) + call mpas_pool_get_config(domain % configs, 'config_use_DMSTracers', config_use_DMSTracers) + call mpas_pool_get_config(domain % configs, 'config_use_MacroMoleculesTracers', config_use_MacroMoleculesTracers) + call mpas_pool_get_config(domain % configs, 'config_use_ecosysTracers_sea_ice_coupling', & + config_use_ecosysTracers_sea_ice_coupling) + call mpas_pool_get_config(domain % configs, 'config_use_DMSTracers_sea_ice_coupling', & + config_use_DMSTracers_sea_ice_coupling) + call mpas_pool_get_config(domain % configs, 'config_use_MacroMoleculesTracers_sea_ice_coupling', & + config_use_MacroMoleculesTracers_sea_ice_coupling) - do i = 1, nCellsSolve - n = n + 1 + n = 0 + block_ptr => domain % blocklist + do while(associated(block_ptr)) + call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) + call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool) - o2x_om(n, index_o2x_So_t) = avgTracersSurfaceValue(index_temperatureSurfaceValue, i) - o2x_om(n, index_o2x_So_s) = avgTracersSurfaceValue(index_salinitySurfaceValue, i) - o2x_om(n, index_o2x_So_u) = avgSurfaceVelocity(index_avgZonalSurfaceVelocity, i) - o2x_om(n, index_o2x_So_v) = avgSurfaceVelocity(index_avgMeridionalSurfaceVelocity, i) + call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) - o2x_om(n, index_o2x_So_ssh) = ssh(i) - o2x_om(n, index_o2x_So_dhdx) = avgSSHGradient(index_avgZonalSSHGradient, i) - o2x_om(n, index_o2x_So_dhdy) = avgSSHGradient(index_avgMeridionalSSHGradient, i) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) - o2x_om(n, index_o2x_Faoo_h2otemp) = avgTotalFreshWaterTemperatureFlux(i) * rho_sw * cp_sw + call mpas_pool_get_dimension(forcingPool, 'index_avgTemperatureSurfaceValue', index_temperatureSurfaceValue) + call mpas_pool_get_dimension(forcingPool, 'index_avgSalinitySurfaceValue', index_salinitySurfaceValue) + call mpas_pool_get_dimension(forcingPool, 'index_avgSurfaceVelocityZonal', index_avgZonalSurfaceVelocity) + call mpas_pool_get_dimension(forcingPool, 'index_avgSurfaceVelocityMeridional', index_avgMeridionalSurfaceVelocity) + call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientZonal', index_avgZonalSSHGradient) + call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientMeridional', index_avgMeridionalSSHGradient) - if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then - o2x_om(n, index_o2x_Foxo_ismw) = avgLandIceFreshwaterFlux(i) - o2x_om(n, index_o2x_Foxo_ismh) = avgLandIceHeatFlux(i) - endif - if (config_remove_ais_river_runoff) then - o2x_om(n, index_o2x_Foxo_rrofl) = avgRemovedRiverRunoffFlux(i) - endif - if (config_remove_ais_ice_runoff) then - o2x_om(n, index_o2x_Foxo_rrofi) = avgRemovedIceRunoffFlux(i) - o2x_om(n, index_o2x_Foxo_rrofih) = avgRemovedIceRunoffHeatFlux(i) - endif + call mpas_pool_get_array(statePool, 'ssh', ssh, 1) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) - if ( frazilIceActive ) then - ! negative when frazil ice can be melted - keepFrazil = .true. - if ( associated(landIceMask) ) then - if ( landIceMask(i) == 1 ) then - keepFrazil = .false. - end if - end if + call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask) + call mpas_pool_get_array(forcingPool, 'avgTracersSurfaceValue', avgTracersSurfaceValue) + call mpas_pool_get_array(forcingPool, 'avgSurfaceVelocity', avgSurfaceVelocity) + call mpas_pool_get_array(forcingPool, 'avgSSHGradient', avgSSHGradient) + call mpas_pool_get_array(forcingPool, 'avgTotalFreshWaterTemperatureFlux', avgTotalFreshWaterTemperatureFlux) + + if ( frazilIceActive ) then + call mpas_pool_get_array(forcingPool, 'seaIceEnergy', seaIceEnergy) + call mpas_pool_get_array(forcingPool, 'frazilSurfacePressure', frazilSurfacePressure) + call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMass, 1) + end if + + ! Cryo fields + if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then + call mpas_pool_get_array(forcingPool, 'avgLandIceFreshwaterFlux', avgLandIceFreshwaterFlux) + call mpas_pool_get_array(forcingPool, 'avgLandIceHeatFlux', avgLandIceHeatFlux) + endif + if (config_remove_ais_river_runoff) then + call mpas_pool_get_array(forcingPool, 'avgRemovedRiverRunoffFlux', avgRemovedRiverRunoffFlux) + endif + if (config_remove_ais_ice_runoff) then + call mpas_pool_get_array(forcingPool, 'avgRemovedIceRunoffFlux', avgRemovedIceRunoffFlux) + call mpas_pool_get_array(forcingPool, 'avgRemovedIceRunoffHeatFlux', avgRemovedIceRunoffHeatFlux) + endif - if ( keepFrazil ) then + ! BGC fields + if (config_use_ecosysTracers) then - ! Calculate energy associated with frazil mass transfer to sea ice if frazil has accumulated - if ( accumulatedFrazilIceMass(i) > 0.0_RKIND ) then + call mpas_pool_get_subpool(forcingPool, 'ecosysAuxiliary', ecosysAuxiliary) + call mpas_pool_get_array(ecosysAuxiliary, 'avgCO2_gas_flux', avgCO2_gas_flux) - seaIceEnergy(i) = accumulatedFrazilIceMass(i) * config_frazil_heat_of_fusion + end if - ! Otherwise calculate the melt potential where avgTracersSurfaceValue represents only the - ! top layer of the ocean - else + if (config_use_ecosysTracers .and. config_use_ecosysTracers_sea_ice_coupling) then + call mpas_pool_get_subpool(forcingPool, 'ecosysSeaIceCoupling', ecosysSeaIceCoupling) - surfaceFreezingTemp = ocn_freezing_temperature(salinity=avgTracersSurfaceValue(index_salinitySurfaceValue, i), & - pressure=0.0_RKIND, inLandIceCavity=.false.) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfacePhytoC', avgOceanSurfacePhytoC) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceDIC', avgOceanSurfaceDIC) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceNO3', avgOceanSurfaceNO3) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceSiO3', avgOceanSurfaceSiO3) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceNH4', avgOceanSurfaceNH4) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceDOCr', avgOceanSurfaceDOCr) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceDOCSemiLabile', avgOceanSurfaceDOCSemiLabile) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceFeParticulate', avgOceanSurfaceFeParticulate) + call mpas_pool_get_array(ecosysSeaIceCoupling, 'avgOceanSurfaceFeDissolved', avgOceanSurfaceFeDissolved) + endif + if (config_use_DMSTracers .and. config_use_DMSTracers_sea_ice_coupling) then + call mpas_pool_get_subpool(forcingPool, 'DMSSeaIceCoupling', DMSSeaIceCoupling) - seaIceEnergy(i) = min(rho_sw*cp_sw*layerThickness(1, i)*( surfaceFreezingTemp + T0_Kelvin & - - avgTracersSurfaceValue(index_temperatureSurfaceValue, i) ), 0.0_RKIND ) + call mpas_pool_get_array(DMSSeaIceCoupling, 'avgOceanSurfaceDMS', avgOceanSurfaceDMS) + call mpas_pool_get_array(DMSSeaIceCoupling, 'avgOceanSurfaceDMSP', avgOceanSurfaceDMSP) + endif + if (config_use_MacroMoleculesTracers .and. config_use_MacroMoleculesTracers_sea_ice_coupling) then + call mpas_pool_get_subpool(forcingPool, 'MacroMoleculesSeaIceCoupling', MacroMoleculesSeaIceCoupling) - end if + call mpas_pool_get_array(MacroMoleculesSeaIceCoupling, 'avgOceanSurfaceDOC', avgOceanSurfaceDOC) + call mpas_pool_get_array(MacroMoleculesSeaIceCoupling, 'avgOceanSurfaceDON', avgOceanSurfaceDON) + endif +! call mpas_pool_get_array(forcingPool, 'CO2Flux', CO2Flux) +! call mpas_pool_get_array(forcingPool, 'DMSFlux', DMSFlux) +! call mpas_pool_get_array(forcingPool, 'surfaceUpwardCO2Flux', surfaceUpwardCO2Flux) - o2x_om(n, index_o2x_Fioo_q) = seaIceEnergy(i) / ocn_cpl_dt - o2x_om(n, index_o2x_Fioo_frazil) = accumulatedFrazilIceMass(i) / ocn_cpl_dt +! replace 'o2x_o % rAttr(' with 'o2x_om(n, ' and ', n)' with ')' + do i = 1, nCellsSolve + n = n + 1 - else + o2x_om(n, index_o2x_So_t) = avgTracersSurfaceValue(index_temperatureSurfaceValue, i) + o2x_om(n, index_o2x_So_s) = avgTracersSurfaceValue(index_salinitySurfaceValue, i) + o2x_om(n, index_o2x_So_u) = avgSurfaceVelocity(index_avgZonalSurfaceVelocity, i) + o2x_om(n, index_o2x_So_v) = avgSurfaceVelocity(index_avgMeridionalSurfaceVelocity, i) - o2x_om(n, index_o2x_Fioo_q) = 0.0_RKIND - o2x_om(n, index_o2x_Fioo_frazil) = 0.0_RKIND - if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then - o2x_om(n, index_o2x_Foxo_q_li) = accumulatedFrazilIceMass(i) * config_frazil_heat_of_fusion / ocn_cpl_dt - o2x_om(n, index_o2x_Foxo_frazil_li) = accumulatedFrazilIceMass(i) / ocn_cpl_dt - endif - end if + o2x_om(n, index_o2x_So_ssh) = ssh(i) + o2x_om(n, index_o2x_So_dhdx) = avgSSHGradient(index_avgZonalSSHGradient, i) + o2x_om(n, index_o2x_So_dhdy) = avgSSHGradient(index_avgMeridionalSSHGradient, i) - ! Reset SeaIce Energy and Accumulated Frazil Ice - seaIceEnergy(i) = 0.0_RKIND - accumulatedFrazilIceMass(i) = 0.0_RKIND - frazilSurfacePressure(i) = 0.0_RKIND - end if + o2x_om(n, index_o2x_Faoo_h2otemp) = avgTotalFreshWaterTemperatureFlux(i) * rho_sw * cp_sw - ! BGC fields - if (config_use_ecosysTracers) then - ! convert from mmolC/m2/s to kg CO2/m2/s - o2x_om(n, index_o2x_Faoo_fco2_ocn) = avgCO2_gas_flux(i)*44.e-6_RKIND - endif - if (config_use_ecosysTracers .and. config_use_ecosysTracers_sea_ice_coupling) then - o2x_om(n, index_o2x_So_algae1) = max(0.0_RKIND,avgOceanSurfacePhytoC(1,i)) - o2x_om(n, index_o2x_So_algae2) = max(0.0_RKIND,avgOceanSurfacePhytoC(2,i)) - o2x_om(n, index_o2x_So_algae3) = max(0.0_RKIND,avgOceanSurfacePhytoC(3,i)) - o2x_om(n, index_o2x_So_dic1) = max(0.0_RKIND,avgOceanSurfaceDIC(i)) - o2x_om(n, index_o2x_So_doc1) = max(0.0_RKIND,avgOceanSurfaceDOCSemiLabile(i)) - o2x_om(n, index_o2x_So_doc2) = max(0.0_RKIND,avgOceanSurfaceDOCSemiLabile(i)) - o2x_om(n, index_o2x_So_doc3) = max(0.0_RKIND,avgOceanSurfaceDOCSemiLabile(i)) - o2x_om(n, index_o2x_So_don1) = 0.0_RKIND - o2x_om(n, index_o2x_So_no3) = max(0.0_RKIND,avgOceanSurfaceNO3(i)) - o2x_om(n, index_o2x_So_sio3) = max(0.0_RKIND,avgOceanSurfaceSiO3(i)) - o2x_om(n, index_o2x_So_nh4) = max(0.0_RKIND,avgOceanSurfaceNH4(i)) - o2x_om(n, index_o2x_So_docr) = max(0.0_RKIND,avgOceanSurfaceDOCr(i)) - o2x_om(n, index_o2x_So_fep1) = max(0.0_RKIND,avgOceanSurfaceFeParticulate(i)) - o2x_om(n, index_o2x_So_fed1) = max(0.0_RKIND,avgOceanSurfaceFeDissolved(i)) - endif - if (config_use_DMSTracers .and. config_use_DMSTracers_sea_ice_coupling) then - o2x_om(n, index_o2x_So_dms) = max(0.0_RKIND,avgOceanSurfaceDMS(i)) - o2x_om(n, index_o2x_So_dmsp) = max(0.0_RKIND,avgOceanSurfaceDMSP(i)) - endif - if (config_use_MacroMoleculesTracers .and. config_use_MacroMoleculesTracers_sea_ice_coupling) then - o2x_om(n, index_o2x_So_doc1) = max(0.0_RKIND,avgOceanSurfaceDOC(1,i)) - o2x_om(n, index_o2x_So_doc2) = max(0.0_RKIND,avgOceanSurfaceDOC(2,i)) - o2x_om(n, index_o2x_So_don1) = max(0.0_RKIND,avgOceanSurfaceDON(i)) - endif + ! Cryo fields + if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then + o2x_om(n, index_o2x_Foxo_ismw) = avgLandIceFreshwaterFlux(i) + o2x_om(n, index_o2x_Foxo_ismh) = avgLandIceHeatFlux(i) + endif + if (config_remove_ais_river_runoff) then + o2x_om(n, index_o2x_Foxo_rrofl) = avgRemovedRiverRunoffFlux(i) + endif + if (config_remove_ais_ice_runoff) then + o2x_om(n, index_o2x_Foxo_rrofi) = avgRemovedIceRunoffFlux(i) + o2x_om(n, index_o2x_Foxo_rrofih) = avgRemovedIceRunoffHeatFlux(i) + endif - if ( trim(config_land_ice_flux_mode) .eq. 'standalone' .or. & - trim(config_land_ice_flux_mode) .eq. 'coupled' ) then - o2x_om(n, index_o2x_So_blt) = landIceBoundaryLayerTracers(indexBLT,i) - o2x_om(n, index_o2x_So_bls) = landIceBoundaryLayerTracers(indexBLS,i) - o2x_om(n, index_o2x_So_htv) = landIceTracerTransferVelocities(indexHeatTrans,i) - o2x_om(n, index_o2x_So_stv) = landIceTracerTransferVelocities(indexSaltTrans,i) - o2x_om(n, index_o2x_So_rhoeff) = 0.0_RKIND - endif - end do + if ( frazilIceActive ) then + ! negative when frazil ice can be melted + keepFrazil = .true. + if ( associated(landIceMask) ) then + if ( landIceMask(i) == 1 ) then + keepFrazil = .false. + end if + end if - block_ptr => block_ptr % next - end do + if ( keepFrazil ) then + + ! Calculate energy associated with frazil mass transfer to sea ice if frazil has accumulated + if ( accumulatedFrazilIceMass(i) > 0.0_RKIND ) then + + seaIceEnergy(i) = accumulatedFrazilIceMass(i) * config_frazil_heat_of_fusion + + ! Otherwise calculate the melt potential where avgTracersSurfaceValue represents only the + ! top layer of the ocean + else + + surfaceFreezingTemp = ocn_freezing_temperature(salinity=avgTracersSurfaceValue(index_salinitySurfaceValue, i), & + pressure=0.0_RKIND, inLandIceCavity=.false.) + seaIceEnergy(i) = min(rho_sw*cp_sw*layerThickness(1, i)*( surfaceFreezingTemp + T0_Kelvin & + - avgTracersSurfaceValue(index_temperatureSurfaceValue, i) ), 0.0_RKIND ) + + end if + + o2x_om(n, index_o2x_Fioo_q) = seaIceEnergy(i) / ocn_cpl_dt + o2x_om(n, index_o2x_Fioo_frazil) = accumulatedFrazilIceMass(i) / ocn_cpl_dt + + else + + o2x_om(n, index_o2x_Fioo_q) = 0.0_RKIND + o2x_om(n, index_o2x_Fioo_frazil) = 0.0_RKIND + if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then + o2x_om(n, index_o2x_Foxo_q_li) = accumulatedFrazilIceMass(i) * config_frazil_heat_of_fusion / ocn_cpl_dt + o2x_om(n, index_o2x_Foxo_frazil_li) = accumulatedFrazilIceMass(i) / ocn_cpl_dt + endif + + end if + + ! Reset SeaIce Energy and Accumulated Frazil Ice + seaIceEnergy(i) = 0.0_RKIND + accumulatedFrazilIceMass(i) = 0.0_RKIND + frazilSurfacePressure(i) = 0.0_RKIND + end if + + ! BGC fields + if (config_use_ecosysTracers .and. index_o2x_Faoo_fco2_ocn /= 0) then + ! convert from mmolC/m2/s to kg CO2/m2/s + o2x_om(n, index_o2x_Faoo_fco2_ocn) = avgCO2_gas_flux(i)*44.e-6_RKIND + endif + if (config_use_ecosysTracers .and. config_use_ecosysTracers_sea_ice_coupling) then + o2x_om(n, index_o2x_So_algae1) = max(0.0_RKIND,avgOceanSurfacePhytoC(1,i)) + o2x_om(n, index_o2x_So_algae2) = max(0.0_RKIND,avgOceanSurfacePhytoC(2,i)) + o2x_om(n, index_o2x_So_algae3) = max(0.0_RKIND,avgOceanSurfacePhytoC(3,i)) + o2x_om(n, index_o2x_So_dic1) = max(0.0_RKIND,avgOceanSurfaceDIC(i)) + o2x_om(n, index_o2x_So_doc1) = max(0.0_RKIND,avgOceanSurfaceDOCSemiLabile(i)) + o2x_om(n, index_o2x_So_doc2) = max(0.0_RKIND,avgOceanSurfaceDOCSemiLabile(i)) + o2x_om(n, index_o2x_So_doc3) = max(0.0_RKIND,avgOceanSurfaceDOCSemiLabile(i)) + o2x_om(n, index_o2x_So_don1) = 0.0_RKIND + o2x_om(n, index_o2x_So_no3) = max(0.0_RKIND,avgOceanSurfaceNO3(i)) + o2x_om(n, index_o2x_So_sio3) = max(0.0_RKIND,avgOceanSurfaceSiO3(i)) + o2x_om(n, index_o2x_So_nh4) = max(0.0_RKIND,avgOceanSurfaceNH4(i)) + o2x_om(n, index_o2x_So_docr) = max(0.0_RKIND,avgOceanSurfaceDOCr(i)) + o2x_om(n, index_o2x_So_fep1) = max(0.0_RKIND,avgOceanSurfaceFeParticulate(i)) + o2x_om(n, index_o2x_So_fed1) = max(0.0_RKIND,avgOceanSurfaceFeDissolved(i)) + endif + if (config_use_DMSTracers .and. config_use_DMSTracers_sea_ice_coupling) then + o2x_om(n, index_o2x_So_dms) = max(0.0_RKIND,avgOceanSurfaceDMS(i)) + o2x_om(n, index_o2x_So_dmsp) = max(0.0_RKIND,avgOceanSurfaceDMSP(i)) + endif + if (config_use_MacroMoleculesTracers .and. config_use_MacroMoleculesTracers_sea_ice_coupling) then + o2x_om(n, index_o2x_So_doc1) = max(0.0_RKIND,avgOceanSurfaceDOC(1,i)) + o2x_om(n, index_o2x_So_doc2) = max(0.0_RKIND,avgOceanSurfaceDOC(2,i)) + o2x_om(n, index_o2x_So_doc3) = max(0.0_RKIND,avgOceanSurfaceDOC(3,i)) + o2x_om(n, index_o2x_So_don1) = 0.0_RKIND + endif +! o2x_om(n, index_o2x_Faoo_fco2_ocn) = CO2Flux(i) +! o2x_om(n, index_o2x_Faoo_fdms_ocn) = DMSFlux(i) +! o2x_om(n, index_o2x_Faoo_fco2_ocn) = surfaceUpwardCO2Flux(i) + +!JW o2x_om(n, index_o2x_So_blt) = landIceBoundaryLayerTemperature(i) +!JW o2x_om(n, index_o2x_So_bls) = landIceBoundaryLayerSalinity(i) +!JW o2x_om(n, index_o2x_So_htv) = landIceHeatTransferVelocity(i) +!JW o2x_om(n, index_o2x_So_stv) = landIceSaltTransferVelocity(i) + + if ( trim(config_land_ice_flux_mode) .eq. 'standalone' .or. & + trim(config_land_ice_flux_mode) .eq. 'coupled' ) then + o2x_om(n, index_o2x_So_blt) = landIceBoundaryLayerTracers(indexBLT,i) + o2x_om(n, index_o2x_So_bls) = landIceBoundaryLayerTracers(indexBLS,i) + o2x_om(n, index_o2x_So_htv) = landIceTracerTransferVelocities(indexHeatTrans,i) + o2x_om(n, index_o2x_So_stv) = landIceTracerTransferVelocities(indexSaltTrans,i) + o2x_om(n, index_o2x_So_rhoeff) = 0.0_RKIND + endif + + !Fyke: test + !write(stderrUnit,*) 'n=',n + !write(stderrUnit,*) 'o2x_om(n, index_o2x_So_blt)=',o2x_om(n, index_o2x_So_blt) + !write(stderrUnit,*) 'o2x_om(n, index_o2x_So_bls)=',o2x_om(n, index_o2x_So_bls) + !write(stderrUnit,*) 'o2x_om(n, index_o2x_So_htv)=',o2x_om(n, index_o2x_So_htv) + !write(stderrUnit,*) 'o2x_om(n, index_o2x_So_stv)=',o2x_om(n, index_o2x_So_stv) + !write(stderrUnit,*) 'o2x_om(n, index_o2x_So_rhoeff)=',o2x_om(n, index_o2x_So_rhoeff) + !o2x_om(n, index_o2x_So_blt) = 0._r8 + !o2x_om(n, index_o2x_So_bls) = 34.5_r8 + !o2x_om(n, index_o2x_So_htv) = 1.e-4_r8 + !o2x_om(n, index_o2x_So_stv) = 3.e-6_r8 + !o2x_om(n, index_o2x_So_rhoeff) = 1000._r8*9.81_r8*918._r8 !lithostatic pressure of 1km of ice + + end do + + block_ptr => block_ptr % next + end do ent_type = 1 ! cells ! set all tags in one method tagname = trim(seq_flds_o2x_fields)//C_NULL_CHAR @@ -4251,9 +4505,12 @@ subroutine ocn_export_moab(EClock) !{{{ wopts = 'PARALLEL=WRITE_PART'//C_NULL_CHAR ierr = iMOAB_WriteMesh(MPOID, outfile, wopts) #endif - end subroutine ocn_export_moab!}}} -#endif +!----------------------------------------------------------------------- +!EOC + + end subroutine ocn_export_moab!}}} +#endif end module ocn_comp_mct !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/driver-moab/main/component_mod.F90 b/driver-moab/main/component_mod.F90 index 3e8ba9042a58..c7a8f51b99ed 100644 --- a/driver-moab/main/component_mod.F90 +++ b/driver-moab/main/component_mod.F90 @@ -471,6 +471,8 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, ! 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: 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 + use seq_comm_mct, only: mb_rof_aream_computed #endif ! ! Arguments @@ -527,6 +529,9 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, dom_s%data%rAttr(km,:) = dom_s%data%rAttr(ka,:) #ifdef HAVE_MOAB + ! TODO should actually compute aream from mesh model + ! we do a lot of unnecessary gymnastics, and very inefficient, because we have a + ! different distribution compared to mct source grid atm tagtype = 1 ! dense, double tagname='aream'//C_NULL_CHAR nloc = mct_avect_lsize(dom_s%data) @@ -542,6 +547,8 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, write(logunit,*) subname,' error in setting the aream tag on atm ' call shr_sys_abort(subname//' ERROR in setting aream tag on atm ') endif + deallocate(gids) + deallocate(data1) ! 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') @@ -597,7 +604,40 @@ subroutine component_init_aream(infodata, rof_c2_ocn, samegrid_ao, samegrid_al, gsmap_s=gsmap_s, av_s=dom_s%data, avfld_s='aream', filefld_s='area_a', & string='rof2ocn ice aream initialization') call t_stopf('CPL:seq_map_readdata-rof2ocn_ice') - + ! this should be more efficient if we just compute aream on coupler side, from actual mesh that we have + ! we need to expose that method in iMOAB, which is local + ! what we do here, we get aream from the domain dom_rx, we just filled it above, with readdata + if(.not. mb_rof_aream_computed) then + + ! we do a lot of unnecessary gymnastics, and very inefficient, because we have a + ! different distribution compared to mct source grid atm + tagtype = 1 ! dense, double + tagname='aream'//C_NULL_CHAR + nloc = mct_avect_lsize(dom_s%data) + allocate(data1(nloc)) + km = mct_aVect_indexRa(dom_s%data, "aream" ) + data1 = dom_s%data%rAttr(km,:) + ent_type = 1 ! element dense double tags + allocate(gids(nloc)) + gids = dom_s%data%iAttr(mct_aVect_indexIA(dom_s%data,"GlobGridNum"),:) + ! ! now set data on the coupler side too + ierr = iMOAB_SetDoubleTagStorageWithGid ( mbrxid, tagname, nloc, ent_type, & + data1, gids) + if (ierr .ne. 0) then + write(logunit,*) subname,' error in setting the aream tag on rof ' + call shr_sys_abort(subname//' ERROR in setting aream tag on rof ') + endif + deallocate(gids) + deallocate(data1) +#ifdef MOABDEBUG + ierr = iMOAB_WriteMesh(mbrxid, trim('recRofWithAream.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 endif end if diff --git a/driver-moab/main/cplcomp_exchange_mod.F90 b/driver-moab/main/cplcomp_exchange_mod.F90 index 76737ad31419..259f6953cff9 100644 --- a/driver-moab/main/cplcomp_exchange_mod.F90 +++ b/driver-moab/main/cplcomp_exchange_mod.F90 @@ -24,7 +24,6 @@ module cplcomp_exchange_mod use seq_comm_mct, only : mhpgid ! iMOAB app id for atm pgx grid, on atm pes use seq_comm_mct, only : atm_pg_active ! flag if PG mesh instanced use seq_comm_mct, only : mlnid , mblxid ! iMOAB app id for land , on land pes and coupler pes - use seq_comm_mct, only : mb_land_mesh ! if true mesh for land use seq_comm_mct, only : mphaid ! iMOAB app id for phys atm; comp atm is 5, phys 5+200 use seq_comm_mct, only : MPSIID, mbixid ! sea-ice on comp pes and on coupler pes use seq_comm_mct, only : mrofid, mbrxid ! iMOAB id of moab rof app on comp pes and on coupler too @@ -1520,16 +1519,12 @@ subroutine cplcomp_moab_Init(infodata,comp) ! we are now on joint pes, compute comm graph between lnd and coupler model typeA = 2 ! point cloud on component PEs, land - if (mb_land_mesh) then - typeA = 3 - endif typeB = 3 ! full mesh on coupler pes, we just read it if (mlnid >= 0) then ierr = iMOAB_GetMeshInfo ( mlnid, nvert, nvise, nbl, nsurf, nvisBC ) comp%mbApCCid = mlnid ! phys atm comp%mbGridType = typeA - 2 ! 0 or 1, pc or cells comp%mblsize = nvert(1) ! vertices - if (mb_land_mesh) comp%mblsize = nvise(1) ! cells endif ierr = iMOAB_ComputeCommGraph( mlnid, mblxid, mpicom_join, mpigrp_old, mpigrp_cplid, & typeA, typeB, id_old, id_join) diff --git a/driver-moab/main/prep_lnd_mod.F90 b/driver-moab/main/prep_lnd_mod.F90 index 7bdb12376a3d..097782f330d6 100644 --- a/driver-moab/main/prep_lnd_mod.F90 +++ b/driver-moab/main/prep_lnd_mod.F90 @@ -18,6 +18,7 @@ module prep_lnd_mod use seq_comm_mct, only: mbrxid ! iMOAB id of moab rof on coupler pes (FV now) use seq_comm_mct, only: mbintxal ! iMOAB id for intx mesh between atm and lnd use seq_comm_mct, only: mbintxrl ! iMOAB id for intx mesh between river and land + use seq_comm_mct, only: mb_rof_aream_computed ! signal use seq_comm_mct, only: mbaxid ! iMOAB id for atm migrated mesh to coupler pes use seq_comm_mct, only: atm_pg_active ! whether the atm uses FV mesh or not ; made true if fv_nphys > 0 @@ -327,6 +328,9 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln fNoBubble, monotonicity, volumetric, fInverseDistanceMap, & noConserve, validate, & trim(dofnameS), trim(dofnameT) ) + + ! signal that the aream for rof has been computed + mb_rof_aream_computed = .true. if (ierr .ne. 0) then write(logunit,*) subname,' error in computing rl weights ' call shr_sys_abort(subname//' ERROR in computing rl weights ') diff --git a/driver-moab/main/prep_rof_mod.F90 b/driver-moab/main/prep_rof_mod.F90 index 83a8e331e916..2ab02f7bdeae 100644 --- a/driver-moab/main/prep_rof_mod.F90 +++ b/driver-moab/main/prep_rof_mod.F90 @@ -30,6 +30,7 @@ module prep_rof_mod use component_type_mod, only: ocn ! used for context for projection towards ocean from rof ! use prep_lnd_mod, only: prep_lnd_get_mapper_Fr2l use map_lnd2rof_irrig_mod, only: map_lnd2rof_irrig + use seq_comm_mct, only: mb_rof_aream_computed ! signal use iso_c_binding #ifdef MOABCOMP @@ -417,6 +418,7 @@ subroutine prep_rof_init(infodata, lnd_c2_rof, atm_c2_rof, ocn_c2_rof) fNoBubble, monotonicity, volumetric, fInverseDistanceMap, & noConserve, validate, & trim(dofnameS), trim(dofnameT) ) + mb_rof_aream_computed = .true. ! signal if (ierr .ne. 0) then write(logunit,*) subname,' error in computing lr weights ' call shr_sys_abort(subname//' ERROR in computing lr weights ') diff --git a/driver-moab/shr/seq_comm_mct.F90 b/driver-moab/shr/seq_comm_mct.F90 index 44d29e91ad2c..9c85df84d2b1 100644 --- a/driver-moab/shr/seq_comm_mct.F90 +++ b/driver-moab/shr/seq_comm_mct.F90 @@ -245,9 +245,9 @@ module seq_comm_mct integer, public :: mbintxar ! iMOAB id for intx mesh between atm and river integer, public :: mbintxlr ! iMOAB id for intx mesh between land and river integer, public :: mbintxrl ! iMOAB id for intx mesh between river and land - logical, public :: mb_land_mesh = .false. ! whether the land uses full FV mesh or not ; made true if domain mesh is read on comp land - + integer, public :: num_moab_exports ! iMOAB id for atm phys grid, on atm pes + logical, public :: mb_rof_aream_computed = .false. ! whether the aream for rof has been set or not !======================================================================= contains