From bba77661a6f0a55c90417b0964ef480c70ae8e70 Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Thu, 22 Jan 2026 10:32:53 +1100 Subject: [PATCH 1/8] some work --- mediator/med_phases_post_rof_mod.F90 | 277 +++++++++++++++++++++++++++ 1 file changed, 277 insertions(+) diff --git a/mediator/med_phases_post_rof_mod.F90 b/mediator/med_phases_post_rof_mod.F90 index 472502f21..de091e05f 100644 --- a/mediator/med_phases_post_rof_mod.F90 +++ b/mediator/med_phases_post_rof_mod.F90 @@ -20,10 +20,13 @@ 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 : field_getdata2d => med_methods_Field_getdata2d use perf_mod , only : t_startf, t_stopf use shr_log_mod , only : shr_log_error + implicit none private @@ -38,8 +41,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,6 +55,8 @@ 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 +74,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,9 +105,19 @@ 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 @@ -128,6 +149,10 @@ subroutine med_phases_post_rof(gcomp, rc) call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) end if + if (spread_rofi) then + call med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) + endif + nullify(is_local%wrap) call ESMF_GridCompGetInternalState(gcomp, is_local, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -161,6 +186,17 @@ 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)) return + 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 + 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 +444,245 @@ 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 + + ! 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), pointer :: areas(:), lats(:) + real(r8), pointer :: rof2ocn_spread(:,:) + real(r8), pointer :: runoff_flux(:) ! temporary 1d pointer + real(r8) :: local_sh(1), global_sh(1) !Antarctic (frozen) runoff + real(r8) :: local_nh(1), global_nh(1) !Greenland (frozen) runoff + + integer :: n, i + + 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)' + !--------------------------------------- + + 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_LogWrite(trim(subname)//": tryign to access gcomp", ESMF_LOGMSG_INFO) + + call ESMF_GridCompGetInternalState(gcomp, is_local, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_LogWrite(trim(subname)//": tryign to get gcomp", ESMF_LOGMSG_INFO) + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! ------------------------------- + ! Create module fields on rof mesh + ! ------------------------------- + + ! if (dbug_flag > dbug_threshold) then + call ESMF_LogWrite(trim(subname)//": tryign to access mesh", ESMF_LOGMSG_INFO) + ! end if + 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 + + do n = 1, size(fields_to_spread_runoff) + + ! assume climatology + ! if (dbug_flag > dbug_threshold) then + call ESMF_LogWrite(trim(subname)//": tryign to create field", ESMF_LOGMSG_INFO) + ! end if + 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 + + call ESMF_LogWrite(trim(subname)//": tryign to read rof2ocn_spread from file", ESMF_LOGMSG_INFO) + call med_io_read(rof2ocn_ice_spread, vm, FBrof_pattern, pre='pattern', rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + do n = 1, size(fields_to_spread_runoff) + + call fldbun_getdata2d(FBrof_pattern, trim(fields_to_spread_runoff(n)), rof2ocn_spread, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + runoff_flux = rof2ocn_spread(:,1) + areas => is_local%wrap%mesh_info(comprof)%areas + lats => is_local%wrap%mesh_info(comprof)%lats + + local_sh(1) = 0.0_r8 + local_nh(1) = 0.0_r8 + do i = 1, size(runoff_flux) + if (lats(i) < 0.0_r8) then + local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) + else + local_nh(1) = local_nh(1) + areas(i) * runoff_flux(i) + end if + end do + + call ESMF_VMAllreduce(vm, senddata=local_sh, recvdata=global_sh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! if (maintask .and. dbug_flag > dbug_threshold) then + if (maintask) then + write(logunit,'(a)') subname//' Correction: '//trim(fields_to_spread_runoff(n)) + write(logunit,'(a,e27.17)') subname//' Correction_sh = ', global_sh(1) + write(logunit,'(a,e27.17)') subname//' Correction_nh = ', global_nh(1) + end if + + ! adjust correction so that it's sums to 1 in each hemisphere + do i = 1, size(runoff_flux) + if (runoff_flux(i) > 0.0000001_r8) then ! is there a better accuracy number to use here ? + if (lats(i) < 0.0_r8) then + runoff_flux(i) = runoff_flux(i) / global_sh(1) + else + runoff_flux(i) = runoff_flux(i) / global_nh(1) + end if + end if + end do + + + local_sh(1) = 0.0_r8 + local_nh(1) = 0.0_r8 + do i = 1, size(runoff_flux) + if (lats(i) < 0.0_r8) then + local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) + else + local_nh(1) = local_nh(1) + areas(i) * runoff_flux(i) + end if + end do + + call ESMF_VMAllreduce(vm, senddata=local_sh, recvdata=global_sh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! if (maintask .and. dbug_flag > dbug_threshold) then + if (maintask) then + write(logunit,'(a)') subname//' Correction: '//trim(fields_to_spread_runoff(n)) + write(logunit,'(a,e27.17)') subname//' Correction_sh = ', global_sh(1) + write(logunit,'(a,e27.17)') subname//' Correction_nh = ', global_nh(1) + end if + + 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 + real(r8), pointer :: runoff_flux(:) ! temporary 1d pointer + real(r8), pointer :: rof2ocn_spread(:,:) + real(r8), pointer :: areas(:), lats(:) + real(r8) :: local_sh(1), global_sh(1) !Antarctic (frozen) runoff + real(r8) :: local_nh(1), global_nh(1) !Greenland (frozen) runoff + real(r8) :: global_sum + ! real(r8) :: multiplier + ! real(r8) :: local_positive_final(1), global_positive_final(1) + ! real(r8) :: local_negative_final(1), global_negative_final(1) + ! real(r8) :: global_sum_final + integer :: n + + 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 + + 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_sh(1) = 0.0_r8 + local_nh(1) = 0.0_r8 + do n = 1, size(runoff_flux) + if (lats(n) < 0.0_r8) then + local_sh(1) = local_sh(1) + areas(n) * runoff_flux(n) + else + local_nh(1) = local_nh(1) + 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_sh, recvdata=global_sh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! if (maintask .and. dbug_flag > dbug_threshold) then + if (maintask) then + write(logunit,'(a)') subname//' Before correction: '//trim(field_name) + write(logunit,'(a,e27.17)') subname//' global_sh = ', global_sh(1) + write(logunit,'(a,e27.17)') subname//' global_nh = ', global_nh(1) + end if + + ! call field_getdata2d(field_rof2ocn_ice_spread, rof2ocn_spread, rc) + ! if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! spear runoff by the saved pattern so that it's sums to 1 in each hemisphere + do n = 1, size(runoff_flux) + if (runoff_flux(n) > 0.0000001_r8) then ! is there a better accuracy number to use here ? + if (lats(n) < 0.0_r8) then + runoff_flux(n) = rof2ocn_spread(n,1) * global_sh(1) + else + runoff_flux(n) = rof2ocn_spread(n,1) * global_nh(1) + end if + end if + end do + + 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 From 2efa4eeff747784d9e5ca2881dba2c984ef2efc2 Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Tue, 24 Feb 2026 10:25:00 +1100 Subject: [PATCH 2/8] working --- mediator/med_io_mod.F90 | 250 ++++++++++++++++++++------- mediator/med_phases_post_rof_mod.F90 | 236 +++++++++++++++---------- 2 files changed, 332 insertions(+), 154 deletions(-) diff --git a/mediator/med_io_mod.F90 b/mediator/med_io_mod.F90 index c86f87c72..f58e90165 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,31 @@ 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 constained 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 + call pio_setframe(pioid, varid, lframe) + 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 +1696,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 +1718,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 +1918,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 +2002,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 +2058,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 de091e05f..e144a70f2 100644 --- a/mediator/med_phases_post_rof_mod.F90 +++ b/mediator/med_phases_post_rof_mod.F90 @@ -4,11 +4,12 @@ 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 @@ -22,7 +23,8 @@ module med_phases_post_rof_mod 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 : field_getdata2d => med_methods_Field_getdata2d + 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 @@ -57,7 +59,7 @@ module med_phases_post_rof_mod 'Forr_rofi_glc'] character(len=9), parameter :: fields_to_spread_runoff(1) = & ['Forr_rofi'] - + character(*) , parameter :: u_FILE_u = & __FILE__ @@ -120,6 +122,7 @@ subroutine med_phases_post_rof_init(gcomp, rc) 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 @@ -139,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)' !--------------------------------------- @@ -149,8 +153,10 @@ subroutine med_phases_post_rof(gcomp, rc) call ESMF_LogWrite(trim(subname)//": called", ESMF_LOGMSG_INFO) end if - if (spread_rofi) then + ! 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) @@ -189,10 +195,16 @@ subroutine med_phases_post_rof(gcomp, rc) 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)) return + 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 @@ -447,6 +459,7 @@ 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 ! input/output variables type(ESMF_GridComp) :: gcomp @@ -457,18 +470,21 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) 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_sh(1), global_sh(1) !Antarctic (frozen) runoff real(r8) :: local_nh(1), global_nh(1) !Greenland (frozen) runoff - - integer :: n, i + 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) @@ -489,101 +505,112 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) ! Create module fields on rof mesh ! ------------------------------- - ! if (dbug_flag > dbug_threshold) then - call ESMF_LogWrite(trim(subname)//": tryign to access mesh", ESMF_LOGMSG_INFO) - ! end if 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 - do n = 1, size(fields_to_spread_runoff) + ! To-do: and support annual/monthly/daily patterns ? + ! for now, assume monthly - ! assume climatology - ! if (dbug_flag > dbug_threshold) then - call ESMF_LogWrite(trim(subname)//": tryign to create field", ESMF_LOGMSG_INFO) - ! end if + ! 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 - + ! call med_methods_FB_Field_diagnose(FBrof_pattern, trim(fields_to_spread_runoff(n)), rc=rc) + ! if (ChkErr(rc,__LINE__,u_FILE_u)) return end do - call ESMF_LogWrite(trim(subname)//": tryign to read rof2ocn_spread from file", ESMF_LOGMSG_INFO) - call med_io_read(rof2ocn_ice_spread, vm, FBrof_pattern, pre='pattern', rc=rc) + ! 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 + ! call med_methods_FB_Field_diagnose(FBrof_pattern, 'Forr_rofi', 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_getdata2d(FBrof_pattern, trim(fields_to_spread_runoff(n)), rof2ocn_spread, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - runoff_flux = rof2ocn_spread(:,1) - areas => is_local%wrap%mesh_info(comprof)%areas - lats => is_local%wrap%mesh_info(comprof)%lats + 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_sh(1) = 0.0_r8 + local_nh(1) = 0.0_r8 + do i = 1, size(rof2ocn_spread, dim = 1) + if (lats(i) < 0.0_r8) then + local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) + else + local_nh(1) = local_nh(1) + areas(i) * runoff_flux(i) + end if + end do - local_sh(1) = 0.0_r8 - local_nh(1) = 0.0_r8 - do i = 1, size(runoff_flux) - if (lats(i) < 0.0_r8) then - local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) - else - local_nh(1) = local_nh(1) + areas(i) * runoff_flux(i) + call ESMF_VMAllreduce(vm, senddata=local_sh, recvdata=global_sh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! if (maintask .and. dbug_flag > dbug_threshold) then + if (maintask) then + write(logunit,'(a)') subname//' Correction: '//trim(fields_to_spread_runoff(n)) + write(logunit,'(a,e27.17)') subname//' Correction_sh = ', global_sh(1) + write(logunit,'(a,e27.17)') subname//' Correction_nh = ', global_nh(1) end if - end do - - call ESMF_VMAllreduce(vm, senddata=local_sh, recvdata=global_sh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! if (maintask .and. dbug_flag > dbug_threshold) then - if (maintask) then - write(logunit,'(a)') subname//' Correction: '//trim(fields_to_spread_runoff(n)) - write(logunit,'(a,e27.17)') subname//' Correction_sh = ', global_sh(1) - write(logunit,'(a,e27.17)') subname//' Correction_nh = ', global_nh(1) - end if - ! adjust correction so that it's sums to 1 in each hemisphere - do i = 1, size(runoff_flux) - if (runoff_flux(i) > 0.0000001_r8) then ! is there a better accuracy number to use here ? + ! adjust correction so that it's sums to 1 in each hemisphere + do i = 1, size(runoff_flux) + if (ieee_is_nan(runoff_flux(i))) then + write(logunit,'(a,e27.17)') subname//' is nan = ', lats(i) + end if + if (runoff_flux(i) .ne. 0) then + if (lats(i) < 0.0_r8) then + runoff_flux(i) = runoff_flux(i) / global_sh(1) + else + runoff_flux(i) = runoff_flux(i) / global_nh(1) + end if + endif + end do + + local_sh(1) = 0.0_r8 + local_nh(1) = 0.0_r8 + do i = 1, size(runoff_flux) if (lats(i) < 0.0_r8) then - runoff_flux(i) = runoff_flux(i) / global_sh(1) + local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) else - runoff_flux(i) = runoff_flux(i) / global_nh(1) + local_nh(1) = local_nh(1) + areas(i) * runoff_flux(i) end if - end if - end do - + end do - local_sh(1) = 0.0_r8 - local_nh(1) = 0.0_r8 - do i = 1, size(runoff_flux) - if (lats(i) < 0.0_r8) then - local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) - else - local_nh(1) = local_nh(1) + areas(i) * runoff_flux(i) + call ESMF_VMAllreduce(vm, senddata=local_sh, recvdata=global_sh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! if (maintask .and. dbug_flag > dbug_threshold) then + if (maintask) then + write(logunit,'(a)') subname//' Correction: '//trim(fields_to_spread_runoff(n)) + write(logunit,'(a,e27.17)') subname//' Correction_sh = ', global_sh(1) + write(logunit,'(a,e27.17)') subname//' Correction_nh = ', global_nh(1) end if - end do - - call ESMF_VMAllreduce(vm, senddata=local_sh, recvdata=global_sh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! if (maintask .and. dbug_flag > dbug_threshold) then - if (maintask) then - write(logunit,'(a)') subname//' Correction: '//trim(fields_to_spread_runoff(n)) - write(logunit,'(a,e27.17)') subname//' Correction_sh = ', global_sh(1) - write(logunit,'(a,e27.17)') subname//' Correction_nh = ', global_nh(1) - end if + rof2ocn_spread(:,month) = runoff_flux + enddo ! month enddo if (dbug_flag > dbug_threshold) then @@ -593,7 +620,6 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) 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. @@ -606,6 +632,8 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, 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(:) @@ -616,7 +644,7 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) ! real(r8) :: local_positive_final(1), global_positive_final(1) ! real(r8) :: local_negative_final(1), global_negative_final(1) ! real(r8) :: global_sum_final - integer :: n + 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)' @@ -629,6 +657,14 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) 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 @@ -664,20 +700,44 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) write(logunit,'(a,e27.17)') subname//' global_nh = ', global_nh(1) end if - ! call field_getdata2d(field_rof2ocn_ice_spread, rof2ocn_spread, rc) - ! if (chkerr(rc,__LINE__,u_FILE_u)) return + !get from fieldbundle + call fldbun_getdata2d(FBrof_pattern, trim(field_name), rof2ocn_spread, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! spear runoff by the saved pattern so that it's sums to 1 in each hemisphere + ! spead runoff by the saved pattern for the model month do n = 1, size(runoff_flux) - if (runoff_flux(n) > 0.0000001_r8) then ! is there a better accuracy number to use here ? - if (lats(n) < 0.0_r8) then - runoff_flux(n) = rof2ocn_spread(n,1) * global_sh(1) - else - runoff_flux(n) = rof2ocn_spread(n,1) * global_nh(1) - end if + if (lats(n) < 0.0_r8) then + runoff_flux(n) = rof2ocn_spread(n,mm) * global_sh(1) + else + runoff_flux(n) = rof2ocn_spread(n,mm) * global_nh(1) end if end do + local_sh(1) = 0.0_r8 + local_nh(1) = 0.0_r8 + do n = 1, size(runoff_flux) + if (lats(n) < 0.0_r8) then + local_sh(1) = local_sh(1) + areas(n) * runoff_flux(n) + else + local_nh(1) = local_nh(1) + 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_sh, recvdata=global_sh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! if (maintask .and. dbug_flag > dbug_threshold) then + if (maintask) then + write(logunit,'(a)') subname//' After correction: '//trim(field_name) + write(logunit,'(a,e27.17)') subname//' global_sh = ', global_sh(1) + write(logunit,'(a,e27.17)') subname//' global_nh = ', global_nh(1) + end if + if (dbug_flag > dbug_threshold) then call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO) end if From df6a11630adb005d1e740087a8d412364955b29c Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Mon, 2 Mar 2026 10:44:25 +1100 Subject: [PATCH 3/8] remove some debugging --- mediator/med_phases_post_rof_mod.F90 | 30 ++++------------------------ 1 file changed, 4 insertions(+), 26 deletions(-) diff --git a/mediator/med_phases_post_rof_mod.F90 b/mediator/med_phases_post_rof_mod.F90 index e144a70f2..776d17f06 100644 --- a/mediator/med_phases_post_rof_mod.F90 +++ b/mediator/med_phases_post_rof_mod.F90 @@ -493,12 +493,10 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) end if nullify(is_local%wrap) - call ESMF_LogWrite(trim(subname)//": tryign to access gcomp", ESMF_LOGMSG_INFO) call ESMF_GridCompGetInternalState(gcomp, is_local, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_LogWrite(trim(subname)//": tryign to get gcomp", ESMF_LOGMSG_INFO) call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! ------------------------------- @@ -521,8 +519,6 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, 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 - ! call med_methods_FB_Field_diagnose(FBrof_pattern, trim(fields_to_spread_runoff(n)), rc=rc) - ! if (ChkErr(rc,__LINE__,u_FILE_u)) return end do ! read spreading from file @@ -531,8 +527,7 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) 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 - ! call med_methods_FB_Field_diagnose(FBrof_pattern, 'Forr_rofi', 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 @@ -565,12 +560,6 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & reduceflag=ESMF_REDUCE_SUM, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! if (maintask .and. dbug_flag > dbug_threshold) then - if (maintask) then - write(logunit,'(a)') subname//' Correction: '//trim(fields_to_spread_runoff(n)) - write(logunit,'(a,e27.17)') subname//' Correction_sh = ', global_sh(1) - write(logunit,'(a,e27.17)') subname//' Correction_nh = ', global_nh(1) - end if ! adjust correction so that it's sums to 1 in each hemisphere do i = 1, size(runoff_flux) @@ -602,14 +591,9 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & reduceflag=ESMF_REDUCE_SUM, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! if (maintask .and. dbug_flag > dbug_threshold) then - if (maintask) then - write(logunit,'(a)') subname//' Correction: '//trim(fields_to_spread_runoff(n)) - write(logunit,'(a,e27.17)') subname//' Correction_sh = ', global_sh(1) - write(logunit,'(a,e27.17)') subname//' Correction_nh = ', global_nh(1) - end if rof2ocn_spread(:,month) = runoff_flux + enddo ! month enddo @@ -640,10 +624,6 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) real(r8) :: local_sh(1), global_sh(1) !Antarctic (frozen) runoff real(r8) :: local_nh(1), global_nh(1) !Greenland (frozen) runoff real(r8) :: global_sum - ! real(r8) :: multiplier - ! real(r8) :: local_positive_final(1), global_positive_final(1) - ! real(r8) :: local_negative_final(1), global_negative_final(1) - ! real(r8) :: global_sum_final integer :: n, mm integer, parameter :: dbug_threshold = 20 ! threshold for writing debug information in this subroutine @@ -693,8 +673,7 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & reduceflag=ESMF_REDUCE_SUM, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! if (maintask .and. dbug_flag > dbug_threshold) then - if (maintask) then + 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_sh(1) write(logunit,'(a,e27.17)') subname//' global_nh = ', global_nh(1) @@ -731,8 +710,7 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & reduceflag=ESMF_REDUCE_SUM, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! if (maintask .and. dbug_flag > dbug_threshold) then - if (maintask) then + if (maintask .and. dbug_flag > dbug_threshold) then write(logunit,'(a)') subname//' After correction: '//trim(field_name) write(logunit,'(a,e27.17)') subname//' global_sh = ', global_sh(1) write(logunit,'(a,e27.17)') subname//' global_nh = ', global_nh(1) From f42d5ed8e228a5011c3983ce9331c9c4f2a3820b Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Mon, 2 Mar 2026 10:59:54 +1100 Subject: [PATCH 4/8] tidyup --- mediator/med_phases_post_rof_mod.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/mediator/med_phases_post_rof_mod.F90 b/mediator/med_phases_post_rof_mod.F90 index 776d17f06..31e73cedc 100644 --- a/mediator/med_phases_post_rof_mod.F90 +++ b/mediator/med_phases_post_rof_mod.F90 @@ -459,7 +459,6 @@ 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 ! input/output variables type(ESMF_GridComp) :: gcomp From 514deb06e5cc3a98174e66a278b22c31fd4fb8ce Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Mon, 2 Mar 2026 11:31:12 +1100 Subject: [PATCH 5/8] fix syntax --- mediator/med_phases_post_rof_mod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/mediator/med_phases_post_rof_mod.F90 b/mediator/med_phases_post_rof_mod.F90 index 31e73cedc..7e7b64ae3 100644 --- a/mediator/med_phases_post_rof_mod.F90 +++ b/mediator/med_phases_post_rof_mod.F90 @@ -459,6 +459,7 @@ 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 From 3fb4f00f7e8ac637153b912fd5ab0f65827c61fe Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Tue, 3 Mar 2026 15:41:41 +1100 Subject: [PATCH 6/8] corrections / neatness --- mediator/med_io_mod.F90 | 6 +++++- mediator/med_phases_post_rof_mod.F90 | 17 ++++++----------- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/mediator/med_io_mod.F90 b/mediator/med_io_mod.F90 index f58e90165..c52758ec0 100644 --- a/mediator/med_io_mod.F90 +++ b/mediator/med_io_mod.F90 @@ -1669,7 +1669,11 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, ungridded_nc, rc) call ESMF_LogWrite(trim(subname)//trim(tmpstr), ESMF_LOGMSG_INFO, rc=rc) fldptr2 = 0.0_r8 else - call pio_setframe(pioid, varid, lframe) + 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) diff --git a/mediator/med_phases_post_rof_mod.F90 b/mediator/med_phases_post_rof_mod.F90 index 7e7b64ae3..4a2885954 100644 --- a/mediator/med_phases_post_rof_mod.F90 +++ b/mediator/med_phases_post_rof_mod.F90 @@ -472,7 +472,7 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) 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 :: rof2ocn_spread(:,:) real(r8), pointer :: runoff_flux(:) ! temporary 1d pointer real(r8) :: local_sh(1), global_sh(1) !Antarctic (frozen) runoff real(r8) :: local_nh(1), global_nh(1) !Greenland (frozen) runoff @@ -546,7 +546,7 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) ! calculate sum of spreading local_sh(1) = 0.0_r8 local_nh(1) = 0.0_r8 - do i = 1, size(rof2ocn_spread, dim = 1) + do i = 1, size(runoff_flux) if (lats(i) < 0.0_r8) then local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) else @@ -563,16 +563,11 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) ! adjust correction so that it's sums to 1 in each hemisphere do i = 1, size(runoff_flux) - if (ieee_is_nan(runoff_flux(i))) then - write(logunit,'(a,e27.17)') subname//' is nan = ', lats(i) + if (lats(i) < 0.0_r8) then + runoff_flux(i) = runoff_flux(i) / global_sh(1) + else + runoff_flux(i) = runoff_flux(i) / global_nh(1) end if - if (runoff_flux(i) .ne. 0) then - if (lats(i) < 0.0_r8) then - runoff_flux(i) = runoff_flux(i) / global_sh(1) - else - runoff_flux(i) = runoff_flux(i) / global_nh(1) - end if - endif end do local_sh(1) = 0.0_r8 From 988c698ed6973c9d3c3cd7a14f497468ba94519f Mon Sep 17 00:00:00 2001 From: anton-seaice Date: Wed, 18 Mar 2026 11:36:51 +1100 Subject: [PATCH 7/8] Some cleanup --- mediator/med_phases_post_rof_mod.F90 | 102 ++++++++++----------------- 1 file changed, 36 insertions(+), 66 deletions(-) diff --git a/mediator/med_phases_post_rof_mod.F90 b/mediator/med_phases_post_rof_mod.F90 index 4a2885954..8fe288d9e 100644 --- a/mediator/med_phases_post_rof_mod.F90 +++ b/mediator/med_phases_post_rof_mod.F90 @@ -12,7 +12,7 @@ module med_phases_post_rof_mod 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 @@ -474,8 +474,7 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) real(r8), pointer :: areas(:), lats(:) real(r8), pointer :: rof2ocn_spread(:,:) real(r8), pointer :: runoff_flux(:) ! temporary 1d pointer - real(r8) :: local_sh(1), global_sh(1) !Antarctic (frozen) runoff - real(r8) :: local_nh(1), global_nh(1) !Greenland (frozen) runoff + real(r8) :: local_sum(2), global_sum(2) ! Antarctic, Greenland (frozen) runoff character(len=CL) :: tempstr integer :: n, i, month @@ -544,49 +543,28 @@ subroutine med_phases_post_rof_init_rof_spread_rofi(gcomp, rc) do month = 1, 12 runoff_flux => rof2ocn_spread(:,month) ! calculate sum of spreading - local_sh(1) = 0.0_r8 - local_nh(1) = 0.0_r8 + local_sum = 0.0_r8 do i = 1, size(runoff_flux) if (lats(i) < 0.0_r8) then - local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) + local_sum(1) = local_sum(1) + areas(i) * runoff_flux(i) else - local_nh(1) = local_nh(1) + areas(i) * runoff_flux(i) + local_sum(2) = local_sum(2) + areas(i) * runoff_flux(i) end if end do - call ESMF_VMAllreduce(vm, senddata=local_sh, recvdata=global_sh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & + 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_sh(1) - else - runoff_flux(i) = runoff_flux(i) / global_nh(1) - end if - end do - - local_sh(1) = 0.0_r8 - local_nh(1) = 0.0_r8 - do i = 1, size(runoff_flux) - if (lats(i) < 0.0_r8) then - local_sh(1) = local_sh(1) + areas(i) * runoff_flux(i) + runoff_flux(i) = runoff_flux(i) / global_sum(1) else - local_nh(1) = local_nh(1) + areas(i) * runoff_flux(i) + runoff_flux(i) = runoff_flux(i) / global_sum(2) end if end do - call ESMF_VMAllreduce(vm, senddata=local_sh, recvdata=global_sh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - rof2ocn_spread(:,month) = runoff_flux enddo ! month @@ -616,9 +594,7 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) real(r8), pointer :: runoff_flux(:) ! temporary 1d pointer real(r8), pointer :: rof2ocn_spread(:,:) real(r8), pointer :: areas(:), lats(:) - real(r8) :: local_sh(1), global_sh(1) !Antarctic (frozen) runoff - real(r8) :: local_nh(1), global_nh(1) !Greenland (frozen) runoff - real(r8) :: global_sum + 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 @@ -650,28 +626,24 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) call fldbun_getdata1d(FBrof_r, trim(field_name), runoff_flux, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - local_sh(1) = 0.0_r8 - local_nh(1) = 0.0_r8 + local_sum = 0.0_r8 do n = 1, size(runoff_flux) if (lats(n) < 0.0_r8) then - local_sh(1) = local_sh(1) + areas(n) * runoff_flux(n) + local_sum(1) = local_sum(1) + areas(n) * runoff_flux(n) else - local_nh(1) = local_nh(1) + areas(n) * runoff_flux(n) + 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_sh, recvdata=global_sh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & + 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_sh(1) - write(logunit,'(a,e27.17)') subname//' global_nh = ', global_nh(1) + 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 @@ -681,34 +653,32 @@ subroutine med_phases_post_rof_spread_rofi(gcomp, field_name, rc) ! 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_sh(1) + runoff_flux(n) = rof2ocn_spread(n,mm) * global_sum(1) else - runoff_flux(n) = rof2ocn_spread(n,mm) * global_nh(1) + runoff_flux(n) = rof2ocn_spread(n,mm) * global_sum(2) end if end do - local_sh(1) = 0.0_r8 - local_nh(1) = 0.0_r8 - do n = 1, size(runoff_flux) - if (lats(n) < 0.0_r8) then - local_sh(1) = local_sh(1) + areas(n) * runoff_flux(n) - else - local_nh(1) = local_nh(1) + 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_sh, recvdata=global_sh, count=1, & - reduceflag=ESMF_REDUCE_SUM, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllreduce(vm, senddata=local_nh, recvdata=global_nh, count=1, & - 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//' After correction: '//trim(field_name) - write(logunit,'(a,e27.17)') subname//' global_sh = ', global_sh(1) - write(logunit,'(a,e27.17)') subname//' global_nh = ', global_nh(1) + + ! 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 From 0153923da97f4f780891f5730540b317fcad9a5a Mon Sep 17 00:00:00 2001 From: Anton Steketee <79179784+anton-seaice@users.noreply.github.com> Date: Wed, 18 Mar 2026 11:37:28 +1100 Subject: [PATCH 8/8] Apply suggestion from @anton-seaice --- mediator/med_io_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mediator/med_io_mod.F90 b/mediator/med_io_mod.F90 index c52758ec0..08fc406c0 100644 --- a/mediator/med_io_mod.F90 +++ b/mediator/med_io_mod.F90 @@ -1659,7 +1659,7 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, ungridded_nc, rc) deallocate(fldptr1_tmp) else if (rank >= 2 .and. lungridded_nc ) then - ! Whole 2d/3d field is constained within one netcdf variable with this name + ! 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)