diff --git a/mediator/med_io_mod.F90 b/mediator/med_io_mod.F90 index c86f87c7..08fc406c 100644 --- a/mediator/med_io_mod.F90 +++ b/mediator/med_io_mod.F90 @@ -7,12 +7,13 @@ module med_io_mod use med_kind_mod , only : CX=>SHR_KIND_CX, CS=>SHR_KIND_CS, CL=>SHR_KIND_CL, I8=>SHR_KIND_I8, R8=>SHR_KIND_R8 use med_kind_mod , only : R4=>SHR_KIND_R4 use med_constants_mod , only : fillvalue => SHR_CONST_SPVAL - use ESMF , only : ESMF_VM, ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_LogFoundError + use ESMF , only : ESMF_VM, ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_LOGMSG_ERROR, ESMF_LogFoundError use ESMF , only : ESMF_SUCCESS, ESMF_LOGERR_PASSTHRU use ESMF , only : ESMF_VMGetCurrent, ESMF_VMGet, ESMF_VMBroadCast, ESMF_Finalize use NUOPC , only : NUOPC_FieldDictionaryGetEntry use NUOPC , only : NUOPC_FieldDictionaryHasEntry use pio , only : file_desc_t, iosystem_desc_t + use pio , only : pio_seterrorhandling, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR use med_internalstate_mod , only : logunit, med_id, maintask use med_constants_mod , only : dbug_flag => med_constants_dbug_flag use med_methods_mod , only : FB_getFieldN => med_methods_FB_getFieldN @@ -20,6 +21,8 @@ module med_io_mod use med_methods_mod , only : FB_getNameN => med_methods_FB_getNameN use med_utils_mod , only : chkerr => med_utils_ChkErr use shr_log_mod , only : shr_log_error + use shr_sys_mod , only : shr_sys_abort + implicit none private @@ -169,6 +172,8 @@ subroutine med_io_init(gcomp, rc) !------------------------------------------------------------------------------- #endif + rc = ESMF_SUCCESS + #ifdef CESMCOUPLED io_subsystem => shr_pio_getiosys(med_id) pio_iotype = shr_pio_getiotype(med_id) @@ -495,10 +500,10 @@ subroutine med_io_wopen(filename, io_file, vm, rc, clobber, file_ind, model_doi_ ! open netcdf file !--------------- - use pio , only : PIO_IOTYPE_PNETCDF, PIO_IOTYPE_NETCDF, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR + use pio , only : PIO_IOTYPE_PNETCDF, PIO_IOTYPE_NETCDF use pio , only : pio_openfile, pio_createfile, PIO_GLOBAL, pio_enddef use pio , only : pio_put_att, pio_redef, pio_get_att - use pio , only : pio_seterrorhandling, pio_file_is_open, pio_clobber, pio_write, pio_noclobber + use pio , only : pio_file_is_open, pio_clobber, pio_write, pio_noclobber ! input/output arguments character(*), intent(in) :: filename @@ -1461,7 +1466,7 @@ subroutine med_io_write_time(io_file, time_val, tbnds, nt, rc) end subroutine med_io_write_time !=============================================================================== - subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) + subroutine med_io_read_FB(filename, vm, FB, pre, frame, ungridded_nc, rc) !--------------- ! Read FB from netcdf file @@ -1472,24 +1477,27 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) use ESMF , only : ESMF_FieldBundleIsCreated, ESMF_FieldBundleGet use ESMF , only : ESMF_FieldGet, ESMF_MeshGet, ESMF_DistGridGet use pio , only : file_desc_T, var_desc_t, io_desc_t, pio_nowrite, pio_openfile - use pio , only : pio_noerr, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR + use pio , only : pio_noerr use pio , only : pio_inq_varid - use pio , only : pio_double, pio_get_att, pio_seterrorhandling, pio_freedecomp, pio_closefile - use pio , only : pio_read_darray, pio_offset_kind, pio_setframe + use pio , only : pio_double, pio_get_att, pio_freedecomp, pio_closefile + use pio , only : pio_read_darray, pio_offset_kind, pio_setframe, pio_strerror ! input/output arguments character(len=*) ,intent(in) :: filename type(ESMF_VM) ,intent(in) :: vm type(ESMF_FieldBundle) ,intent(in) :: FB ! data to be read character(len=*) ,optional ,intent(in) :: pre ! prefix to variable name + logical ,optional ,intent(in) :: ungridded_nc + ! if true : ungridded dim in fields is dimension in netcdf file, + ! if false: ungridded_dim is provided in seperate variables with the index appended integer(kind=PIO_OFFSET_KIND) ,optional ,intent(in) :: frame integer ,intent(out) :: rc ! local variables type(ESMF_Field) :: lfield - integer :: rcode + integer :: rcode, ierr integer :: nf - integer :: k,n,l + integer :: k,n,l,m type(file_desc_t) :: pioid type(var_desc_t) :: varid type(io_desc_t) :: iodesc @@ -1503,12 +1511,13 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) character(CL) :: tmpstr character(len=16) :: cnumber integer(kind=Pio_Offset_Kind) :: lframe - integer :: ungriddedUBound(1) ! currently the size must equal 1 for rank 2 fieldds - integer :: gridToFieldMap(1) ! currently the size must equal 1 for rank 2 fieldds + integer :: ungriddedUBound(1) ! currently the size must equal 1 for rank 2 fields + integer :: gridToFieldMap(1) ! currently the size must equal 1 for rank 2 fields + logical :: lungridded_nc character(*),parameter :: subName = '(med_io_read_FB) ' !------------------------------------------------------------------------------- rc = ESMF_Success - call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) + call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return lpre = ' ' @@ -1520,11 +1529,16 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) else lframe = 1 endif + if (present(ungridded_nc)) then + lungridded_nc = ungridded_nc + else + lungridded_nc = .false. + endif if (.not. ESMF_FieldBundleIsCreated(FB,rc=rc)) then call ESMF_LogWrite(trim(subname)//" FB "//trim(lpre)//" not created", ESMF_LOGMSG_INFO) if (chkerr(rc,__LINE__,u_FILE_u)) return if (dbug_flag > 5) then - call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) + call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return endif return @@ -1533,10 +1547,10 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) call ESMF_FieldBundleGet(FB, fieldCount=nf, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return write(tmpstr,*) subname//' field count = '//trim(lpre),nf - call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO) + call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (nf < 1) then - call ESMF_LogWrite(trim(subname)//" FB "//trim(lpre)//" empty", ESMF_LOGMSG_INFO) + call ESMF_LogWrite(trim(subname)//" FB "//trim(lpre)//" empty", ESMF_LOGMSG_INFO, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (dbug_flag > 5) then call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) @@ -1547,7 +1561,7 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) if (med_io_file_exists(vm, trim(filename))) then rcode = pio_openfile(io_subsystem, pioid, pio_iotype, trim(filename),pio_nowrite) - call ESMF_LogWrite(trim(subname)//' open file '//trim(filename), ESMF_LOGMSG_INFO) + call ESMF_LogWrite(trim(subname)//' open file '//trim(filename), ESMF_LOGMSG_INFO, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return else call shr_log_error(trim(subname)//' ERROR: file invalid '//trim(filename), & @@ -1568,17 +1582,20 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) call ESMF_FieldBundleGet(FB, itemc, field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldGet(lfield, rank=rank, rc=rc) + if (lungridded_nc .and. rank == 1) then + call shr_sys_abort(trim(subname)//' ERROR: ungridded_nc = true but field only contains one dimension', rc=rc) + endif if (chkerr(rc,__LINE__,u_FILE_u)) return - if (rank == 2) then + if (rank==2 .and. .not. lungridded_nc) then name1 = trim(lpre)//'_'//trim(itemc)//'1' - else if (rank == 1) then + else name1 = trim(lpre)//'_'//trim(itemc) end if - call med_io_read_init_iodesc(FB, name1, pioid, iodesc, rc) + call med_io_read_init_iodesc(lfield, trim(name1), pioid, iodesc, ungridded_nc=lungridded_nc, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return end if - call ESMF_LogWrite(trim(subname)//' reading field '//trim(itemc), ESMF_LOGMSG_INFO) + call ESMF_LogWrite(trim(subname)//' reading field '//trim(itemc), ESMF_LOGMSG_INFO, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! Get pointer to field bundle field @@ -1587,7 +1604,7 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) fldptr1=fldptr1, fldptr2=fldptr2, rank=rank, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - if (rank == 2) then + if (rank == 2 .and. .not. lungridded_nc ) then ! Determine the size of the ungridded dimension and the ! index where the undistributed dimension is located @@ -1611,18 +1628,25 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) rcode = pio_inq_varid(pioid, trim(name1), varid) if (rcode == pio_noerr) then - call ESMF_LogWrite(trim(subname)//' read field '//trim(name1), ESMF_LOGMSG_INFO) + call ESMF_LogWrite(trim(subname)//' read field '//trim(name1), ESMF_LOGMSG_INFO, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call pio_setframe(pioid, varid, lframe) call pio_read_darray(pioid, varid, iodesc, fldptr1_tmp, rcode) - rcode = pio_get_att(pioid, varid, "_FillValue", lfillvalue) if (rcode /= pio_noerr) then - lfillvalue = fillvalue + call ESMF_LogWrite(trim(subname)//' failed to read variable '//trim(name1), ESMF_LOGMSG_INFO, rc=rc) + ierr = pio_strerror(rcode, tmpstr) + call ESMF_LogWrite(trim(subname)//trim(tmpstr), ESMF_LOGMSG_ERROR, rc=rc) + else + rcode = pio_get_att(pioid, varid, "_FillValue", lfillvalue) + if (rcode /= pio_noerr) then + lfillvalue = fillvalue + endif + where (fldptr1_tmp == lfillvalue) fldptr1_tmp = 0.0_r8 endif - do l = 1,size(fldptr1_tmp) - if (fldptr1_tmp(l) == lfillvalue) fldptr1_tmp(l) = 0.0_r8 - enddo else + call ESMF_LogWrite(trim(subname)//' failed to read variable '//trim(name1), ESMF_LOGMSG_INFO, rc=rc) + ierr = pio_strerror(rcode, tmpstr) + call ESMF_LogWrite(trim(subname)//trim(tmpstr), ESMF_LOGMSG_ERROR, rc=rc) fldptr1_tmp = 0.0_r8 endif if (gridToFieldMap(1) == 1) then @@ -1634,6 +1658,35 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) deallocate(fldptr1_tmp) + else if (rank >= 2 .and. lungridded_nc ) then + ! Whole 2d/3d field is contained within one netcdf variable with this name + name1 = trim(lpre)//'_'//trim(itemc) + + rcode = pio_inq_varid(pioid, trim(name1), varid) + if (rcode /= pio_noerr) then + call ESMF_LogWrite(trim(subname)//' failed to read variable '//trim(name1), ESMF_LOGMSG_INFO, rc=rc) + ierr = pio_strerror(rcode, tmpstr) + call ESMF_LogWrite(trim(subname)//trim(tmpstr), ESMF_LOGMSG_INFO, rc=rc) + fldptr2 = 0.0_r8 + else + if (present(frame)) then + call pio_seterrorhandling(pioid,PIO_INTERNAL_ERROR) + call pio_setframe(pioid, varid, lframe) + call pio_seterrorhandling(pioid,PIO_BCAST_ERROR) + endif + call pio_read_darray(pioid, varid, iodesc, fldptr2, rcode) + if (rcode /= pio_noerr) then + call ESMF_LogWrite(trim(subname)//' failed to read variable '//trim(name1), ESMF_LOGMSG_INFO, rc=rc) + ierr = pio_strerror(rcode, tmpstr) + call ESMF_LogWrite(trim(subname)//trim(tmpstr), ESMF_LOGMSG_ERROR, rc=rc) + else + rcode = pio_get_att(pioid, varid, "_FillValue", lfillvalue) + if (rcode /= pio_noerr) then + lfillvalue = fillvalue + endif + where (fldptr2 == lfillvalue) fldptr2 = 0.0_r8 + endif + endif else if (rank == 1) then name1 = trim(lpre)//'_'//trim(itemc) @@ -1647,13 +1700,14 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) if (rcode /= pio_noerr) then lfillvalue = fillvalue endif - do n = 1,size(fldptr1) - if (fldptr1(n) == lfillvalue) fldptr1(n) = 0.0_r8 - enddo + where (fldptr1 == lfillvalue) fldptr1 = 0.0_r8 else fldptr1 = 0.0_r8 endif - end if + else + write(tmpstr,*) rank + call shr_log_error(trim(subname)//': rank='//trim(tmpstr)//' of field '//trim(itemc)//' not supported', rc=rc) + end if enddo ! end of loop over fields call pio_seterrorhandling(pioid,PIO_INTERNAL_ERROR) @@ -1668,101 +1722,169 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc) end subroutine med_io_read_FB !=============================================================================== - subroutine med_io_read_init_iodesc(FB, name1, pioid, iodesc, rc) + subroutine med_io_read_init_iodesc(field, name1, pioid, iodesc, ungridded_nc, rc) use ESMF , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS use ESMF , only : ESMF_FieldBundleIsCreated, ESMF_FieldBundle, ESMF_Mesh, ESMF_DistGrid use ESMF , only : ESMF_FieldBundleGet, ESMF_FieldGet, ESMF_MeshGet, ESMF_DistGridGet use ESMF , only : ESMF_Field, ESMF_FieldGet, ESMF_AttributeGet use pio , only : file_desc_T, var_desc_t, io_desc_t, pio_nowrite, pio_openfile - use pio , only : pio_noerr, pio_inq_varndims + use pio , only : pio_noerr, pio_strerror, pio_inq_varndims use pio , only : pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, pio_inq_vardimid - use pio , only : pio_double, pio_seterrorhandling, pio_initdecomp + use pio , only : pio_double, pio_initdecomp ! input/output variables - type(ESMF_FieldBundle) , intent(in) :: FB + type(ESMF_Field) , intent(in) :: field character(len=*) , intent(in) :: name1 - type(file_desc_t) , intent(in) :: pioid + type(file_desc_t) , intent(inout) :: pioid type(io_desc_t) , intent(inout) :: iodesc + logical ,optional ,intent(in) :: ungridded_nc integer , intent(out) :: rc ! local variables - type(ESMF_Field) :: field type(ESMF_Mesh) :: mesh type(ESMF_Distgrid) :: distgrid - integer :: rcode - integer :: ns,ng + integer :: rcode, ierr + integer :: ns,i integer :: ndims integer, pointer :: dimid(:) type(var_desc_t) :: varid - integer :: lnx,lny - integer, pointer :: minIndexPTile(:,:) + type(integer), allocatable :: gdims(:) + integer :: lnx,lny,lni integer, pointer :: maxIndexPTile(:,:) integer :: dimCount, tileCount integer, pointer :: Dof(:) + logical :: lungridded_nc character(CL) :: tmpstr character(*),parameter :: subName = '(med_io_read_init_iodesc) ' !------------------------------------------------------------------------------- rc = ESMF_SUCCESS + if (present(ungridded_nc)) then + lungridded_nc = ungridded_nc + else + lungridded_nc = .false. + endif rcode = pio_inq_varid(pioid, trim(name1), varid) if (rcode == pio_noerr) then rcode = pio_inq_varndims(pioid, varid, ndims) + if (rcode /= pio_noerr) then + ierr = pio_strerror(rcode, tmpstr) + call shr_log_error(trim(subname)//' ERROR: '//trim(tmpstr), & + line=__LINE__, file=u_FILE_u, rc=rc) + return + endif + write(tmpstr,*) trim(subname),' ndims = ',ndims call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO) + if (ndims>3 .or. ndims <1) then + write(tmpstr,*) ndims + call shr_sys_abort(trim(subname)//' ERROR: ndims = '//trim(tmpstr)//' is not supported. ', rc=rc) + endif + allocate(dimid(ndims)) rcode = pio_inq_vardimid(pioid, varid, dimid(1:ndims)) + if (rcode /= pio_noerr) then + ierr = pio_strerror(rcode, tmpstr) + call shr_log_error(trim(subname)//' ERROR: '//trim(tmpstr), & + line=__LINE__, file=u_FILE_u, rc=rc) + return + endif + rcode = pio_inq_dimlen(pioid, dimid(1), lnx) + if (rcode /= pio_noerr) then + ierr = pio_strerror(rcode, tmpstr) + call shr_sys_abort(trim(subname)//' ERROR: '//trim(tmpstr), rc=rc) + endif write(tmpstr,*) trim(subname),' lnx = ',lnx call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO) - if (ndims>=2) then - rcode = pio_inq_dimlen(pioid, dimid(2), lny) + if ( (ndims==2 .and. .not. lungridded_nc) .or. ndims>2 ) then !2nd dimension is gridded + rcode = pio_inq_dimlen(pioid, dimid(2), lny) + if (rcode /= pio_noerr) then + ierr = pio_strerror(rcode, tmpstr) + call shr_log_error(trim(subname)//' ERROR: '//trim(tmpstr), & + line=__LINE__, file=u_FILE_u, rc=rc) + return + endif + write(tmpstr,*) trim(subname),' lny = ',lny + call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO) else - lny = 1 - end if + lny = 1 !default + endif + if (lungridded_nc) then !2nd/3rd dimension is ungridded + rcode = pio_inq_dimlen(pioid, dimid(ndims), lni) + if (rcode /= pio_noerr) then + ierr = pio_strerror(rcode, tmpstr) + call shr_log_error(trim(subname)//' ERROR: '//trim(tmpstr), & + line=__LINE__, file=u_FILE_u, rc=rc) + return + endif + write(tmpstr,*) trim(subname),' lni = ',lni + call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO) + else + lni = 1 !default + endif deallocate(dimid) - write(tmpstr,*) trim(subname),' lny = ',lny - call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO) - ng = lnx * lny - call FB_getFieldN(FB, 1, field, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + ! get the number of dimensions in the grid call ESMF_FieldGet(field, mesh=mesh, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_MeshGet(mesh, elementDistgrid=distgrid, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_DistGridGet(distgrid, dimCount=dimCount, tileCount=tileCount, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - - allocate(minIndexPTile(dimCount, tileCount), maxIndexPTile(dimCount, tileCount)) - call ESMF_DistGridGet(distgrid, minIndexPTile=minIndexPTile, & - maxIndexPTile=maxIndexPTile, rc=rc) + ! check that grid and data size match + allocate(maxIndexPTile(dimCount, tileCount)) + call ESMF_DistGridGet(distgrid, maxIndexPTile=maxIndexPTile, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - - if (ng > maxval(maxIndexPTile)) then + if ((lnx * lny) > maxval(maxIndexPTile)) then write(tmpstr,*) subname,' WARNING: dimensions do not match', lnx, lny, maxval(maxIndexPTile) call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO) ! This should not be an error for say CTSM which does not send a global grid endif + deallocate(maxIndexPTile) + ! create PIO decomposition + ! get number of elements in local grid call ESMF_DistGridGet(distgrid, localDE=0, elementCount=ns, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - allocate(dof(ns)) + if ( lungridded_nc ) then + allocate(dof(ns*lni)) + else + allocate(dof(ns)) + endif + ! get 1d indeces of local grid on global grid (dof) call ESMF_DistGridGet(distgrid, localDE=0, seqIndexList=dof, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + if ( lungridded_nc ) then + ! dof is index for i=1, cycle around all i (as i is ungridded) + do i=2,lni + dof((i-1)*ns+1:i*ns) = dof(1:ns) + (i-1)*lnx*lny + enddo + endif write(tmpstr,*) subname,' dof = ',ns,size(dof),dof(1),dof(ns) !,minval(dof),maxval(dof) call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO) - call pio_initdecomp(io_subsystem, pio_double, (/lnx,lny/), dof, iodesc) - deallocate(dof) + allocate(gdims(ndims)) + if (ndims==1) gdims = (/lnx/) + if (ndims==2 .and. .not. lungridded_nc) gdims = (/lnx, lny/) + if (ndims==2 .and. lungridded_nc) gdims = (/lnx, lni/) + if (ndims==3 .and. lni > 1) gdims = (/lnx,lny,lni/) + if (ndims==3 .and. lni == 1) gdims = (/lnx,lny/) - deallocate(minIndexPTile, maxIndexPTile) + call pio_seterrorhandling(pioid,PIO_INTERNAL_ERROR) + call pio_initdecomp(io_subsystem, pio_double, gdims, dof, iodesc) + call pio_seterrorhandling(pioid, PIO_BCAST_ERROR) + + deallocate(dof, gdims) else - call shr_log_error(trim(subname)//' ERROR: '//trim(name1)//' is not present, aborting ', rc=rc) - return + ierr = pio_strerror(rcode, tmpstr) + call shr_sys_abort(trim(subname)//' ERROR: '//trim(name1)//' is not present, aborting. '//trim(tmpstr), rc=rc) end if ! end if rcode check end subroutine med_io_read_init_iodesc @@ -1800,7 +1922,7 @@ subroutine med_io_read_int1d(filename, vm, idata, dname, rc) ! Read 1d integer array from netcdf file !--------------- - use pio , only : var_desc_t, file_desc_t, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR, pio_seterrorhandling + use pio , only : var_desc_t use pio , only : pio_get_var, pio_inq_varid, pio_get_att, pio_openfile use pio , only : pio_nowrite, pio_openfile, pio_global use pio , only : pio_closefile @@ -1884,8 +2006,8 @@ subroutine med_io_read_r81d(filename, vm, rdata, dname, rc) ! Read 1d double array from netcdf file !--------------- - use pio , only : file_desc_t, var_desc_t, pio_openfile, pio_closefile, pio_seterrorhandling - use pio , only : PIO_BCAST_ERROR, PIO_INTERNAL_ERROR, pio_inq_varid, pio_get_var + use pio , only : var_desc_t, pio_openfile, pio_closefile + use pio , only : pio_inq_varid, pio_get_var use pio , only : pio_nowrite, pio_openfile, pio_global, pio_get_att ! input/output arguments @@ -1940,7 +2062,7 @@ subroutine med_io_read_char(filename, vm, rdata, dname, rc) ! Read char string from netcdf file !--------------- - use pio , only : file_desc_t, var_desc_t, pio_seterrorhandling, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR + use pio , only : var_desc_t use pio , only : pio_closefile, pio_inq_varid, pio_get_var use pio , only : pio_openfile, pio_global, pio_get_att, pio_nowrite diff --git a/mediator/med_phases_post_rof_mod.F90 b/mediator/med_phases_post_rof_mod.F90 index 472502f2..8fe288d9 100644 --- a/mediator/med_phases_post_rof_mod.F90 +++ b/mediator/med_phases_post_rof_mod.F90 @@ -4,14 +4,15 @@ module med_phases_post_rof_mod use NUOPC_Mediator , only : NUOPC_MediatorGet use NUOPC , only : NUOPC_CompAttributeGet - use ESMF , only : ESMF_Clock, ESMF_ClockIsCreated - use ESMF , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS + use ESMF , only : ESMF_Clock, ESMF_Time, ESMF_ClockIsCreated + use ESMF , only : ESMF_ClockGet, ESMF_TimeGet, ESMF_ClockIsCreated + use ESMF , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS, ESMF_LOGMSG_ERROR use ESMF , only : ESMF_GridComp, ESMF_GridCompGet use ESMF , only : ESMF_Mesh, ESMF_MESHLOC_ELEMENT, ESMF_TYPEKIND_R8 - use ESMF , only : ESMF_Field, ESMF_FieldCreate + use ESMF , only : ESMF_Field, ESMF_FieldCreate, ESMF_FieldGet use ESMF , only : ESMF_FieldBundle, ESMF_FieldBundleCreate use ESMF , only : ESMF_FieldBundleGet, ESMF_FieldBundleAdd - use ESMF , only : ESMF_VM, ESMF_VMAllreduce, ESMF_REDUCE_SUM + use ESMF , only : ESMF_VM, ESMF_VMAllreduce, ESMF_VMReduce, ESMF_REDUCE_SUM use med_kind_mod , only : CX=>SHR_KIND_CX, CS=>SHR_KIND_CS, CL=>SHR_KIND_CL, R8=>SHR_KIND_R8 use med_internalstate_mod , only : complnd, compocn, compice, comprof use med_internalstate_mod , only : InternalState, maintask, logunit @@ -20,10 +21,14 @@ module med_phases_post_rof_mod use med_phases_history_mod, only : med_phases_history_write_comp use med_map_mod , only : med_map_field_packed use med_methods_mod , only : fldbun_getdata1d => med_methods_FB_getdata1d + use med_methods_mod , only : fldbun_getdata2d => med_methods_FB_getdata2d use med_methods_mod , only : fldbun_getmesh => med_methods_FB_getmesh + use med_methods_mod , only : fldbun_getFldPtr => med_methods_FB_getFldPtr + use med_methods_mod , only : FB_fldchk => med_methods_FB_FldChk use perf_mod , only : t_startf, t_stopf use shr_log_mod , only : shr_log_error + implicit none private @@ -38,8 +43,13 @@ module med_phases_post_rof_mod integer :: num_rof_fields character(len=CS), allocatable :: rof_field_names(:) + ! A local FieldBundle to store the distribution(s) pattern for runoff. + type(ESMF_FieldBundle) :: FBrof_pattern + logical :: remove_negative_runoff_lnd logical :: remove_negative_runoff_glc + logical :: spread_rofi + character(len=CL) :: rof2ocn_ice_spread character(len=9), parameter :: fields_to_remove_negative_runoff_lnd(2) = & ['Forr_rofl', & @@ -47,7 +57,9 @@ module med_phases_post_rof_mod character(len=13), parameter :: fields_to_remove_negative_runoff_glc(2) = & ['Forr_rofl_glc', & 'Forr_rofi_glc'] - + character(len=9), parameter :: fields_to_spread_runoff(1) = & + ['Forr_rofi'] + character(*) , parameter :: u_FILE_u = & __FILE__ @@ -64,6 +76,7 @@ subroutine med_phases_post_rof_init(gcomp, rc) ! local variables character(CL) :: cvalue logical :: isPresent, isSet + integer :: n character(len=*), parameter :: subname='(med_phases_post_rof_init)' !--------------------------------------- @@ -94,11 +107,22 @@ subroutine med_phases_post_rof_init(gcomp, rc) remove_negative_runoff_glc = .false. end if + call NUOPC_CompAttributeGet(gcomp, name='rof2ocn_ice_spread', value=rof2ocn_ice_spread, isPresent=isPresent, isSet=isSet, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + spread_rofi = .true. + else + spread_rofi = .false. + end if + if (maintask) then write(logunit,'(a,l7)') trim(subname)//' remove_negative_runoff_lnd = ', remove_negative_runoff_lnd write(logunit,'(a,l7)') trim(subname)//' remove_negative_runoff_glc = ', remove_negative_runoff_glc + write(logunit,'(a,l7)') trim(subname)//' spread_rofi = ', spread_rofi + if (spread_rofi) write(logunit,'(a)') trim(subname)//' rof2ocn_ice_spread = '//trim(rof2ocn_ice_spread) end if + if (dbug_flag > 20) then call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) end if @@ -118,6 +142,7 @@ subroutine med_phases_post_rof(gcomp, rc) real(r8), pointer :: data_copy(:) integer :: n logical :: exists + logical :: first_time = .true. character(len=*), parameter :: subname='(med_phases_post_rof)' !--------------------------------------- @@ -128,6 +153,12 @@ subroutine med_phases_post_rof(gcomp, rc) call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) end if + ! unclear why this can't be in med_phases_post_rof_init, possibly pio not initialised + if (spread_rofi .and. first_time) then + call med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) + first_time=.false. + endif + nullify(is_local%wrap) call ESMF_GridCompGetInternalState(gcomp, is_local, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -161,6 +192,23 @@ subroutine med_phases_post_rof(gcomp, rc) end do end if + if (spread_rofi) then + do n = 1, size(fields_to_spread_runoff) + call ESMF_FieldBundleGet(FBrof_r, fieldName=trim(fields_to_spread_runoff(n)), isPresent=exists, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) then + call shr_log_error(string=trim(subname)//" Error checking field: "//trim(fields_to_spread_runoff(n)), line=__LINE__,file=u_FILE_u, rc=rc) + return + end if + if (exists) then + call med_phases_post_rof_spread_rofi(gcomp, fields_to_spread_runoff(n), rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call shr_log_error(string=trim(subname)//" Runoff field to spread: "//trim(fields_to_spread_runoff(n))//" does not exist", line=__LINE__,file=u_FILE_u, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + end do + end if + ! map rof to lnd if (is_local%wrap%med_coupling_active(comprof,complnd)) then call t_startf('MED:'//trim(subname)//' map_rof2lnd') @@ -408,4 +456,236 @@ subroutine med_phases_post_rof_remove_negative_runoff(gcomp, field_name, rc) end subroutine med_phases_post_rof_remove_negative_runoff + subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) + !--------------------------------------------------------------- + use med_io_mod , only : med_io_read + use, intrinsic :: ieee_arithmetic, only : ieee_is_nan + + ! input/output variables + type(ESMF_GridComp) :: gcomp + integer, intent(out) :: rc + + + type(ESMF_Mesh) :: mesh_l + type(InternalState) :: is_local + type(ESMF_VM) :: vm + type(ESMF_field) :: field_l ! climatology, 12 months + real(r8) :: glob_area_inv + real(r8), pointer :: areas(:), lats(:) + real(r8), pointer :: rof2ocn_spread(:,:) + real(r8), pointer :: runoff_flux(:) ! temporary 1d pointer + real(r8) :: local_sum(2), global_sum(2) ! Antarctic, Greenland (frozen) runoff + character(len=CL) :: tempstr + integer :: n, i, month + + integer, parameter :: dbug_threshold = 0 ! threshold for writing debug information in this subroutine + character(len=*), parameter :: subname='(med_phases_post_rof_mod: med_phases_post_rof_init_rof_spread_rofi)' + !--------------------------------------- + + ! to do - make component configurable (could be comprof or compatm) + + rc = ESMF_SUCCESS + + call t_startf('MED:'//subname) + if (dbug_flag > dbug_threshold) then + call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) + end if + + nullify(is_local%wrap) + + call ESMF_GridCompGetInternalState(gcomp, is_local, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! ------------------------------- + ! Create module fields on rof mesh + ! ------------------------------- + + call fldbun_getmesh(is_local%wrap%FBImp(comprof,comprof), mesh_l, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + FBrof_pattern = ESMF_FieldBundleCreate(name='FBrof_pattern', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! To-do: and support annual/monthly/daily patterns ? + ! for now, assume monthly + + ! create field for each field_to_spread + do n = 1, size(fields_to_spread_runoff) + field_l = ESMF_FieldCreate(mesh_l, ESMF_TYPEKIND_R8, name=trim(fields_to_spread_runoff(n)), meshloc=ESMF_MESHLOC_ELEMENT, & + ungriddedLBound=(/1/), ungriddedUBound=(/12/), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldBundleAdd(FBrof_pattern, (/field_l/), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end do + + ! read spreading from file + if (dbug_flag > dbug_threshold) then + call ESMF_LogWrite(trim(subname)//": tryign to read rof2ocn_spread from file", ESMF_LOGMSG_INFO) + endif + call med_io_read(rof2ocn_ice_spread, vm, FBrof_pattern, pre='pattern', ungridded_nc=.true., rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + areas => is_local%wrap%mesh_info(comprof)%areas + lats => is_local%wrap%mesh_info(comprof)%lats + + ! normalise each month in each hemisphere of the spreading pattern + do n = 1, size(fields_to_spread_runoff) + + call fldbun_getFldPtr(FBrof_pattern, trim(fields_to_spread_runoff(n)), & + fldptr2=rof2ocn_spread, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) then + call ESMF_LogWrite(trim(subname)//": rof2ocn_spread not retrieved", ESMF_LOGMSG_ERROR) + return + endif + + do month = 1, 12 + runoff_flux => rof2ocn_spread(:,month) + ! calculate sum of spreading + local_sum = 0.0_r8 + do i = 1, size(runoff_flux) + if (lats(i) < 0.0_r8) then + local_sum(1) = local_sum(1) + areas(i) * runoff_flux(i) + else + local_sum(2) = local_sum(2) + areas(i) * runoff_flux(i) + end if + end do + + call ESMF_VMAllreduce(vm, senddata=local_sum, recvdata=global_sum, count=2, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! adjust correction so that it's sums to 1 in each hemisphere + do i = 1, size(runoff_flux) + if (lats(i) < 0.0_r8) then + runoff_flux(i) = runoff_flux(i) / global_sum(1) + else + runoff_flux(i) = runoff_flux(i) / global_sum(2) + end if + end do + + rof2ocn_spread(:,month) = runoff_flux + + enddo ! month + enddo + + if (dbug_flag > dbug_threshold) then + call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) + end if + call t_stopf('MED:'//subname) + + end subroutine med_phases_post_rof_init_rof_spread_rofi + + subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) + !--------------------------------------------------------------- + ! For one runoff field, spread runoff according to the pattern prescribed in spread_rofi_weights. + + ! input/output variables + type(ESMF_GridComp) :: gcomp + character(len=*), intent(in) :: field_name ! name of runoff flux field to process + integer, intent(out) :: rc + + ! local variables + type(InternalState) :: is_local + type(ESMF_VM) :: vm + type(ESMF_Clock) :: clock + type(ESMF_Time) :: currTime + real(r8), pointer :: runoff_flux(:) ! temporary 1d pointer + real(r8), pointer :: rof2ocn_spread(:,:) + real(r8), pointer :: areas(:), lats(:) + real(r8) :: local_sum(2), global_sum(2) !Antarctic,Greenland (frozen) runoff + integer :: n, mm + + integer, parameter :: dbug_threshold = 20 ! threshold for writing debug information in this subroutine + character(len=*), parameter :: subname='(med_phases_post_rof_mod: med_phases_post_rof_spread_rofi)' + !--------------------------------------- + + rc = ESMF_SUCCESS + + call t_startf('MED:'//subname) + if (dbug_flag > dbug_threshold) then + call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) + end if + + !get the month of year + call ESMF_GridCompGet(gcomp, clock=clock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_ClockGet(clock, currTime=currTime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_TimeGet(currTime, mm=mm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + nullify(is_local%wrap) + call ESMF_GridCompGetInternalState(gcomp, is_local, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + areas => is_local%wrap%mesh_info(comprof)%areas + lats => is_local%wrap%mesh_info(comprof)%lats + + call fldbun_getdata1d(FBrof_r, trim(field_name), runoff_flux, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + local_sum = 0.0_r8 + do n = 1, size(runoff_flux) + if (lats(n) < 0.0_r8) then + local_sum(1) = local_sum(1) + areas(n) * runoff_flux(n) + else + local_sum(2) = local_sum(2) + areas(n) * runoff_flux(n) + end if + end do + + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_sum, recvdata=global_sum, count=2, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (maintask .and. dbug_flag > dbug_threshold) then + write(logunit,'(a)') subname//' Before correction: '//trim(field_name) + write(logunit,'(a,e27.17)') subname//' global_sh = ', global_sum(1) + write(logunit,'(a,e27.17)') subname//' global_nh = ', global_sum(2) + end if + + !get from fieldbundle + call fldbun_getdata2d(FBrof_pattern, trim(field_name), rof2ocn_spread, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! spead runoff by the saved pattern for the model month + do n = 1, size(runoff_flux) + if (lats(n) < 0.0_r8) then + runoff_flux(n) = rof2ocn_spread(n,mm) * global_sum(1) + else + runoff_flux(n) = rof2ocn_spread(n,mm) * global_sum(2) + end if + end do + + if (maintask .and. dbug_flag > dbug_threshold) then + + ! calculate the new global sum (after correction), should be equal to 0 + local_sum = 0.0_r8 + do n = 1, size(runoff_flux) + if (lats(n) < 0.0_r8) then + local_sum(1) = local_sum(1) + areas(n) * runoff_flux(n) + else + local_sum(2) = local_sum(2) + areas(n) * runoff_flux(n) + end if + end do + + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMReduce(vm, senddata=local_sum, recvdata=global_sum, count=2, & + reduceflag=ESMF_REDUCE_SUM, rootPet = 0, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + write(logunit,'(a)') subname//' Before correction: '//trim(field_name) + write(logunit,'(a,e27.17)') subname//' global_sh = ', global_sum(1) + write(logunit,'(a,e27.17)') subname//' global_nh = ', global_sum(2) + end if + + if (dbug_flag > dbug_threshold) then + call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) + end if + call t_stopf('MED:'//subname) + + end subroutine med_phases_post_rof_spread_rofi + end module med_phases_post_rof_mod