diff --git a/src/ICE/iceupdate.F90 b/src/ICE/iceupdate.F90 index a973312f04364781465403618f175cf52c6170ed..785a68a6979508a76e0cae3cabe397e3b8ac16ba 100644 --- a/src/ICE/iceupdate.F90 +++ b/src/ICE/iceupdate.F90 @@ -31,7 +31,7 @@ MODULE iceupdate USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE lbclnk ! lateral boundary conditions (or mpp links) USE timing ! Timing - + USE oce IMPLICIT NONE PRIVATE @@ -446,12 +446,15 @@ CONTAINS CALL iom_get( numrir, jpdom_auto, 'snwice_mass_b', snwice_mass_b ) ELSE ! start from rest IF(lwp) WRITE(numout,*) ' ==>> previous run without snow-ice mass output then set it' - snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) ) + snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) & + & + rhow * (vt_ip(:,:) + vt_il(:,:)) ) snwice_mass_b(:,:) = snwice_mass(:,:) ENDIF ELSE !* Start from rest +!JC: I think this is useless with what is now done in ice_istate IF(lwp) WRITE(numout,*) ' ==>> start from rest: set the snow-ice mass' - snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) ) + snwice_mass (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) & + & + rhow * (vt_ip(:,:) + vt_il(:,:)) ) snwice_mass_b(:,:) = snwice_mass(:,:) ENDIF ! diff --git a/src/NST/agrif_all_update.F90 b/src/NST/agrif_all_update.F90 index 8fc617eba676216116c165d70709a00c11f1b8e4..58d595d6fd7404275cb677c18a2a80b158d4dad0 100644 --- a/src/NST/agrif_all_update.F90 +++ b/src/NST/agrif_all_update.F90 @@ -42,6 +42,11 @@ CONTAINS ! IF (lwp.AND.lk_agrif_debug) Write(*,*) ' --> START AGRIF UPDATE from grid Number',Agrif_Fixed() ! + ! Update computational domain mask once: + IF (lk_agrif_fstep) THEN + CALL Agrif_Update_Variable(tmask_id,locupdate=(/ nn_shift_bar,-2/), procname = update_tmask_agrif) + ENDIF + ! CALL Agrif_Update_ssh() ! Update sea level ! IF (.NOT.ln_linssh) CALL Agrif_Update_vvl() ! Update scale factors diff --git a/src/NST/agrif_ice_update.F90 b/src/NST/agrif_ice_update.F90 index d142e9ab77275f1233f07d188420e8696e2465e2..4557ef513ef4fa52c4f0e9bef05b97039c5ac8a4 100644 --- a/src/NST/agrif_ice_update.F90 +++ b/src/NST/agrif_ice_update.F90 @@ -24,6 +24,7 @@ MODULE agrif_ice_update USE sbc_oce USE agrif_oce USE ice + USE sbc_ice , ONLY : snwice_mass USE agrif_ice USE phycst , ONLY: rt0 @@ -101,21 +102,21 @@ CONTAINS IF( before ) THEN jm = 1 DO jl = 1, jpl - ptab(i1:i2,j1:j2,jm ) = a_i (i1:i2,j1:j2,jl) - ptab(i1:i2,j1:j2,jm+1) = v_i (i1:i2,j1:j2,jl) - ptab(i1:i2,j1:j2,jm+2) = v_s (i1:i2,j1:j2,jl) - ptab(i1:i2,j1:j2,jm+3) = sv_i(i1:i2,j1:j2,jl) - ptab(i1:i2,j1:j2,jm+4) = oa_i(i1:i2,j1:j2,jl) - ptab(i1:i2,j1:j2,jm+5) = a_ip(i1:i2,j1:j2,jl) - ptab(i1:i2,j1:j2,jm+6) = v_ip(i1:i2,j1:j2,jl) - ptab(i1:i2,j1:j2,jm+7) = v_il(i1:i2,j1:j2,jl) - ptab(i1:i2,j1:j2,jm+8) = t_su(i1:i2,j1:j2,jl) + ptab(i1:i2,j1:j2,jm ) = a_i (i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) + ptab(i1:i2,j1:j2,jm+1) = v_i (i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) + ptab(i1:i2,j1:j2,jm+2) = v_s (i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) + ptab(i1:i2,j1:j2,jm+3) = sv_i(i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) + ptab(i1:i2,j1:j2,jm+4) = oa_i(i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) + ptab(i1:i2,j1:j2,jm+5) = a_ip(i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) + ptab(i1:i2,j1:j2,jm+6) = v_ip(i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) + ptab(i1:i2,j1:j2,jm+7) = v_il(i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) + ptab(i1:i2,j1:j2,jm+8) = t_su(i1:i2,j1:j2,jl) * e1e2t_frac(i1:i2,j1:j2) jm = jm + 9 DO jk = 1, nlay_s - ptab(i1:i2,j1:j2,jm) = e_s(i1:i2,j1:j2,jk,jl) ; jm = jm + 1 + ptab(i1:i2,j1:j2,jm) = e_s(i1:i2,j1:j2,jk,jl) * e1e2t_frac(i1:i2,j1:j2) ; jm = jm + 1 END DO DO jk = 1, nlay_i - ptab(i1:i2,j1:j2,jm) = e_i(i1:i2,j1:j2,jk,jl) ; jm = jm + 1 + ptab(i1:i2,j1:j2,jm) = e_i(i1:i2,j1:j2,jk,jl) * e1e2t_frac(i1:i2,j1:j2) ; jm = jm + 1 END DO END DO ! @@ -164,6 +165,17 @@ CONTAINS DO jl = 1, jpl WHERE( tmask(i1:i2,j1:j2,1) == 0._wp ) t_su(i1:i2,j1:j2,jl) = rt0 ! to avoid a division by 0 in sbcblk.F90 END DO + + ! new mass per unit area + vt_s (i1:i2,j1:j2) = SUM( v_s(i1:i2,j1:j2,:), dim=3 ) + vt_i (i1:i2,j1:j2) = SUM( v_i(i1:i2,j1:j2,:), dim=3 ) + vt_ip(i1:i2,j1:j2) = SUM( v_ip(i1:i2,j1:j2,:), dim=3 ) + vt_il(i1:i2,j1:j2) = SUM( v_il(i1:i2,j1:j2,:), dim=3 ) + + snwice_mass(i1:i2,j1:j2) = tmask(i1:i2,j1:j2,1) * ( rhos * vt_s(i1:i2,j1:j2) & + & + rhoi * vt_i(i1:i2,j1:j2) & + & + rhow * (vt_ip(i1:i2,j1:j2) & + & + vt_il(i1:i2,j1:j2)) ) ENDIF ! diff --git a/src/NST/agrif_oce.F90 b/src/NST/agrif_oce.F90 index 01d8034632dc960c82e8c04363a562fdc5837bd7..691e418458b9c7635bc5806bdf9c04e8d38eaeda 100644 --- a/src/NST/agrif_oce.F90 +++ b/src/NST/agrif_oce.F90 @@ -80,7 +80,7 @@ MODULE agrif_oce INTEGER, PUBLIC :: unb_interp_id, vnb_interp_id, ub2b_interp_id, vb2b_interp_id INTEGER, PUBLIC :: ub2b_update_id, vb2b_update_id, unb_update_id, vnb_update_id INTEGER, PUBLIC :: ub2b_cor_id, vb2b_cor_id, sshn_id - INTEGER, PUBLIC :: sshn_frc_id + INTEGER, PUBLIC :: sshn_frc_id, tmask_id INTEGER, PUBLIC :: e3t_id, e3u_id, e3v_id, e3f_id INTEGER, PUBLIC :: r3t_id, r3u_id, r3v_id, r3f_id INTEGER, PUBLIC :: scales_t_id diff --git a/src/NST/agrif_oce_interp.F90 b/src/NST/agrif_oce_interp.F90 index dc459b771739b8fd3283e7883a4fbd905fcd7614..d4fd3787597faea5238cc3265e5a119b825026b5 100644 --- a/src/NST/agrif_oce_interp.F90 +++ b/src/NST/agrif_oce_interp.F90 @@ -53,6 +53,7 @@ MODULE agrif_oce_interp PUBLIC interp_e1e2t_frac, interp_e2u_frac, interp_e1v_frac PUBLIC agrif_istate_oce, agrif_istate_ssh ! called by icestate.F90 and domvvl.F90 PUBLIC agrif_check_bat + PUBlIC interp_tmask_agrif INTEGER :: bdy_tinterp = 0 @@ -1118,6 +1119,26 @@ CONTAINS END SUBROUTINE interpsshn_frc + SUBROUTINE interp_tmask_agrif( ptab, i1, i2, j1, j2, before ) + !!---------------------------------------------------------------------- + !! *** ROUTINE interp_tmask_agrif *** + !! + !! set tmask_agrif = 0 over ghost points + !! + !!---------------------------------------------------------------------- + INTEGER , INTENT(in ) :: i1, i2, j1, j2 + REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab + LOGICAL , INTENT(in ) :: before + ! + !!---------------------------------------------------------------------- + ! + IF(.NOT.before) THEN + tmask_agrif(i1:i2,j1:j2) = 0._wp + ENDIF + ! + END SUBROUTINE interp_tmask_agrif + + SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE interpun *** diff --git a/src/NST/agrif_oce_update.F90 b/src/NST/agrif_oce_update.F90 index 0f538408182df07010a6dfc9620eac779c7fd4d9..378acd29c7d727da9cdcc9937fbfb9d7235095e2 100644 --- a/src/NST/agrif_oce_update.F90 +++ b/src/NST/agrif_oce_update.F90 @@ -36,6 +36,7 @@ MODULE agrif_oce_update PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh PUBLIC Agrif_Check_parent_bat + PUBLIC update_tmask_agrif !! * Substitutions # include "domzgr_substitute.h90" @@ -1903,6 +1904,24 @@ CONTAINS ENDIF ! END SUBROUTINE check_parent_e3v0 + + + SUBROUTINE update_tmask_agrif( tabres, i1, i2, j1, j2, before ) + !!---------------------------------------------------------------------- + !! *** ROUTINE updatetmsk *** + !!---------------------------------------------------------------------- + INTEGER , INTENT(in ) :: i1, i2, j1, j2 + REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres + LOGICAL , INTENT(in ) :: before + !! + !!---------------------------------------------------------------------- + ! + IF( .NOT.before ) THEN + tmask_agrif(i1:i2,j1:j2) = 0._wp + ENDIF + ! + END SUBROUTINE update_tmask_agrif + #else !!---------------------------------------------------------------------- !! Empty module no AGRIF zoom diff --git a/src/NST/agrif_user.F90 b/src/NST/agrif_user.F90 index 321c5b91d2ab8ceabf81420df48d8f5151b819fd..4af2c2ecfd9c206880b2815c9857a013c4faa647 100644 --- a/src/NST/agrif_user.F90 +++ b/src/NST/agrif_user.F90 @@ -99,6 +99,7 @@ CALL agrif_declare_variable((/1,1 /),(/ind2-1,ind3-1 /),(/'x','y' /),(/1,1 /),(/jpi,jpj /), r3f_id) #endif CALL agrif_declare_variable((/2,2,0 /),(/ind2 ,ind3 ,0 /),(/'x','y','N' /),(/1,1,1 /),(/jpi,jpj,jpk /),e3t0_interp_id) + CALL agrif_declare_variable((/2,2 /),(/ind2 ,ind3 /),(/'x','y' /),(/1,1 /),(/jpi,jpj /), tmask_id) CALL agrif_declare_variable((/2,2 /),(/ind2 ,ind3 /),(/'x','y' /),(/1,1 /),(/jpi,jpj /), mbkt_id) CALL agrif_declare_variable((/2,2 /),(/ind2 ,ind3 /),(/'x','y' /),(/1,1 /),(/jpi,jpj /), ht0_id) CALL agrif_declare_variable((/2,2 /),(/ind2 ,ind3 /),(/'x','y' /),(/1,1 /),(/jpi,jpj /), e1e2t_frac_id) @@ -127,6 +128,7 @@ #endif CALL Agrif_Set_bcinterp(e3t0_interp_id,interp =AGRIF_linear ) CALL Agrif_Set_interp (e3t0_interp_id,interp =AGRIF_linear ) + CALL Agrif_Set_bcinterp( tmask_id,interp =AGRIF_constant) CALL Agrif_Set_bcinterp( mbkt_id,interp =AGRIF_constant) CALL Agrif_Set_interp ( mbkt_id,interp =AGRIF_constant) CALL Agrif_Set_bcinterp( ht0_id,interp =AGRIF_constant) @@ -168,6 +170,7 @@ ! extend the interpolation zone by 1 more point than necessary: ! RB check here CALL Agrif_Set_bc( e3t0_interp_id, (/-nn_sponge_len*imaxrho-2,ind1-1/) ) + CALL Agrif_Set_bc( tmask_id, (/-imaxrho*nn_shift_bar,ind1-1/) ) CALL Agrif_Set_bc( mbkt_id, (/-nn_sponge_len*imaxrho-2,ind1-1/) ) CALL Agrif_Set_bc( ht0_id, (/-nn_sponge_len*imaxrho-2,ind1-1/) ) CALL Agrif_Set_bc( e1e2t_frac_id, (/-nn_sponge_len*imaxrho-2,ind1-1/) ) @@ -198,6 +201,7 @@ CALL Agrif_Set_Updatetype( r3f_id,update = Agrif_Update_Copy ) #endif #endif + CALL Agrif_Set_Updatetype( tmask_id,update = AGRIF_Update_Average) CALL Agrif_Set_ExternalMapping(nemo_mapping) ! @@ -242,6 +246,8 @@ ht0_parent( :,:) = 0._wp mbkt_parent(:,:) = 0 ! + ! Build tmask_agrif such that it is zero outside barotropic dynamical interface: + CALL Agrif_Bc_variable(tmask_id ,calledweight=1.,procname=interp_tmask_agrif) ! CALL Agrif_Bc_variable(ht0_id ,calledweight=1.,procname=interpht0 ) ! CALL Agrif_Bc_variable(mbkt_id,calledweight=1.,procname=interpmbkt) CALL Agrif_Init_Variable(ht0_id, procname=interpht0 ) @@ -410,29 +416,17 @@ hbdy(:,:) = 0._wp ssh(:,:,Krhs_a) = 0._wp - IF ( ln_dynspg_ts ) THEN - Agrif_UseSpecialValue = ln_spc_dyn - use_sign_north = .TRUE. - sign_north = -1. - CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) ! must be called before unb_id to define ubdy - CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) ! must be called before vnb_id to define vbdy - CALL Agrif_Bc_variable( unb_interp_id,calledweight=1.,procname=interpunb ) - CALL Agrif_Bc_variable( vnb_interp_id,calledweight=1.,procname=interpvnb ) - use_sign_north = .FALSE. - ubdy(:,:) = 0._wp - vbdy(:,:) = 0._wp - ELSEIF ( ln_dynspg_EXP ) THEN - Agrif_UseSpecialValue = ln_spc_dyn - use_sign_north = .TRUE. - sign_north = -1. - ubdy(:,:) = 0._wp - vbdy(:,:) = 0._wp - CALL Agrif_Bc_variable( unb_interp_id,calledweight=1.,procname=interpunb ) - CALL Agrif_Bc_variable( vnb_interp_id,calledweight=1.,procname=interpvnb ) - use_sign_north = .FALSE. - ubdy(:,:) = 0._wp - vbdy(:,:) = 0._wp - ENDIF + Agrif_UseSpecialValue = ln_spc_dyn + use_sign_north = .TRUE. + sign_north = -1. + ubdy(:,:) = 0._wp + vbdy(:,:) = 0._wp + CALL Agrif_Bc_variable( unb_interp_id,calledweight=1.,procname=interpunb ) + CALL Agrif_Bc_variable( vnb_interp_id,calledweight=1.,procname=interpvnb ) + use_sign_north = .FALSE. + ubdy(:,:) = 0._wp + vbdy(:,:) = 0._wp + Agrif_UseSpecialValue = .FALSE. l_vremap = .FALSE. diff --git a/src/OCE/DOM/dom_oce.F90 b/src/OCE/DOM/dom_oce.F90 index 7d73a76d14d0446a62b183b4c1d1beaf83bc553a..d61ac9bb059de44ef249a9b3c9318b473d7b6a7e 100644 --- a/src/OCE/DOM/dom_oce.F90 +++ b/src/OCE/DOM/dom_oce.F90 @@ -199,7 +199,8 @@ MODULE dom_oce REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask !: land/ocean mask at T-, U-, V-, W- and F-pts REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wumask, wvmask !: land/ocean mask at WU- and WV-pts REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: fe3mask !: land/ocean mask at F-pts - REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: land/ocean mask at F-pts + REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: agrif mask defining barotropic update + REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_agrif !: agrif mask at T-points excluding ghosts and updated areas !!---------------------------------------------------------------------- !! calendar variables @@ -342,7 +343,8 @@ CONTAINS ! #if defined key_agrif ii = ii+1 - ALLOCATE( tmask_upd(jpi,jpj) , umask_upd(jpi,jpj), vmask_upd(jpi,jpj) , STAT=ierr(ii) ) + ALLOCATE( tmask_upd(jpi,jpj) , umask_upd(jpi,jpj), vmask_upd(jpi,jpj), & + & tmask_agrif(jpi,jpj), STAT=ierr(ii) ) #endif ! dom_oce_alloc = MAXVAL(ierr) diff --git a/src/OCE/DOM/dommsk.F90 b/src/OCE/DOM/dommsk.F90 index 904fec9b03f78748a6aa4e12da5f2420d1f7428d..362b11e3e48d4bfad0279a14b9b72c8847554aeb 100644 --- a/src/OCE/DOM/dommsk.F90 +++ b/src/OCE/DOM/dommsk.F90 @@ -158,6 +158,15 @@ CONTAINS ! In case of a coarsened grid, account her for possibly aditionnal ! masked points; these have been read in the mesh file and stored in mbku, mbkv, mbkf DO_2D( 0, 0, 0, 0 ) + ! Ugly patch to accomodate batropic case (jpk=2) + ! but with possible masked faces in the coarsening case + ! A better way would be to keep track of mbku/v/f=0 until here + ! but these have been assigned a minimum of 1. TBC + IF (jpk>2) THEN + IF (mbku(ji,jj)==1) umask(ji,jj,:) = 0._wp + IF (mbkv(ji,jj)==1) vmask(ji,jj,:) = 0._wp + IF (mbkf(ji,jj)==1) fmask(ji,jj,:) = 0._wp + ENDIF IF ( MAXVAL(umask(ji,jj,:))/=0._wp ) umask(ji,jj,mbku(ji,jj)+1:jpk) = 0._wp IF ( MAXVAL(vmask(ji,jj,:))/=0._wp ) vmask(ji,jj,mbkv(ji,jj)+1:jpk) = 0._wp IF ( MAXVAL(fmask(ji,jj,:))/=0._wp ) fmask(ji,jj,mbkf(ji,jj)+1:jpk) = 0._wp @@ -225,6 +234,10 @@ CONTAINS tmask_upd(:,:) = 0._wp umask_upd(:,:) = 0._wp vmask_upd(:,:) = 0._wp + ! + ! Reset mask defining actual computationnal domain + ! e.g. excluding ghosts and updated cells. + tmask_agrif(:,:) = 1._wp #endif ! END SUBROUTINE dom_msk diff --git a/src/OCE/DYN/dynspg_ts.F90 b/src/OCE/DYN/dynspg_ts.F90 index f8856aedf4774b83c2f8cc215cf73454827f0540..c41e6f167c04cc45edd7ce6c2ea69b5be256b7e3 100644 --- a/src/OCE/DYN/dynspg_ts.F90 +++ b/src/OCE/DYN/dynspg_ts.F90 @@ -1016,6 +1016,12 @@ CONTAINS zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) END_2D ! +#if defined key_agrif + ! Discard points that are not stepped by 2d mode: + zcu(Nis0:Nie0,Njs0:Nje0) = zcu(Nis0:Nie0,Njs0:Nje0) & + & * (1._wp - tmask_upd(Nis0:Nie0,Njs0:Nje0)) +#endif + ! zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) ) CALL mpp_max( 'dynspg_ts', zcmax ) diff --git a/src/OCE/SBC/sbcfwb.F90 b/src/OCE/SBC/sbcfwb.F90 index 58a17f028f4a7cc404859b8c7c4a66f98cd04406..91754a5373b68d3a282971f39c84d61fe79118ee 100644 --- a/src/OCE/SBC/sbcfwb.F90 +++ b/src/OCE/SBC/sbcfwb.F90 @@ -21,6 +21,9 @@ MODULE sbcfwb USE phycst ! physical constants USE sbcrnf ! ocean runoffs USE sbcssr ! Sea-Surface damping terms +#if defined key_agrif + USE agrif_oce , ONLY: Kmm_a +#endif ! USE in_out_manager ! I/O manager USE iom ! IOM @@ -37,8 +40,13 @@ MODULE sbcfwb REAL(wp) :: rn_fwb0 ! initial freshwater adjustment flux [kg/m2/s] (nn_fwb = 2 only) REAL(wp) :: a_fwb ! annual domain averaged freshwater budget from the previous year REAL(wp) :: a_fwb_b ! annual domain averaged freshwater budget from the year before or at initial state - REAL(wp) :: a_fwb_ini ! initial domain averaged freshwater budget REAL(wp) :: area ! global mean ocean surface (interior domain) + REAL(wp) :: emp_corr ! current, globally constant, emp correction +#if defined key_agrif +!$AGRIF_DO_NOT_TREAT + REAL(wp) :: agrif_sum ! global sum used for agrif +!$AGRIF_END_DO_NOT_TREAT +#endif !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) @@ -67,8 +75,9 @@ CONTAINS INTEGER, INTENT( in ) :: Kmm ! ocean time level index ! INTEGER :: ios, inum, ikty ! local integers - REAL(wp) :: z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp ! local scalars - REAL(wp) :: zsurf_neg, zsurf_pos, zsurf_tospread, zcoef ! - - + INTEGER :: ji, jj, istart, iend, jstart, jend + REAL(wp) :: z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp ! local scalars + REAL(wp) :: zsurf_neg, zsurf_pos, zsurf_tospread ! - - REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmsk_neg, ztmsk_pos, z_wgt ! 2D workspaces REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmsk_tospread, zerp_cor ! - - REAL(wp) ,DIMENSION(1) :: z_fwfprv @@ -99,8 +108,20 @@ CONTAINS ! IF( kn_fwb == 3 .AND. nn_sssr /= 2 ) CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' ) IF( kn_fwb == 3 .AND. ln_isfcav ) CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' ) +#if defined key_agrif + IF((kn_fwb == 3).AND.(Agrif_maxlevel()/=0)) CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 not yet implemented with AGRIF zooms ' ) +#endif ! - area = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask(:,:,1)) ! interior global domain surface + IF ( Agrif_Root() ) THEN +#if defined key_agrif + agrif_sum = 0.0_wp ! set work scalar to 0 + CALL Agrif_step_child(glob_sum_area_agrif) ! integrate over grids + area = agrif_sum +#else + area = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask(:,:,1)) ! interior global domain surface +#endif + IF (lwp) WRITE(numout,*) 'Domain area (km**2):', area/1000._wp/1000._wp + ENDIF ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes ! and in case of no melt, it can generate HSSW. ! @@ -116,77 +137,173 @@ CONTAINS ! CASE ( 1 ) !== global mean fwf set to zero ==! ! +#if ! defined key_agrif IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN - y_fwfnow(1) = local_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) ) - CALL mpp_delay_sum( 'sbcfwb', 'fwb', y_fwfnow(:), z_fwfprv(:), kt == nitend - nn_fsbc + 1 ) - z_fwfprv(1) = z_fwfprv(1) / area - zcoef = z_fwfprv(1) * rcp - emp(:,:) = emp(:,:) - z_fwfprv(1) * tmask(:,:,1) - qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction - ! outputs - IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', zcoef * sst_m(:,:) * tmask(:,:,1) ) - IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', z_fwfprv(1) * tmask(:,:,1) ) + z_fwfprv(1) = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) & + & - fwfisf_cav(:,:) - fwfisf_par(:,:) & + & - snwice_fmass(:,:) ) ) / area + emp_corr = -z_fwfprv(1) + emp(:,:) = emp(:,:) + emp_corr * tmask(:,:,1) + qns(:,:) = qns(:,:) - emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction + ENDIF +#else + IF ( Agrif_Root() ) THEN + IF ( Agrif_maxlevel()==0 ) THEN + ! No child grid, correct "now" fluxes (i.e. as in the "no agrif" case) + IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN + emp_corr = -glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) & + & - fwfisf_cav(:,:) - fwfisf_par(:,:) & + & - snwice_fmass(:,:) ) ) / area + ENDIF + ELSE + ! + ! Volume is here corrected according to the budget computed in the past, e.g. between + ! the last two consecutive calls to the surface module. Hence, the volume is allowed to drift slightly during + ! the current time step. + ! + IF( kt == nit000 ) THEN ! initialisation + ! ! + IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb', ldstop = .FALSE. ) > 0 & ! read from restart file + & .AND. iom_varid( numror, 'emp_corr', ldstop = .FALSE. ) > 0 ) THEN + IF(lwp) WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file' + CALL iom_get( numror, 'a_fwb' , a_fwb ) + CALL iom_get( numror, 'emp_corr' , emp_corr ) + ! + ELSE + emp_corr = 0._wp + a_fwb = 999._wp + a_fwb_b = 999._wp + END IF + ! + IF(lwp) WRITE(numout,*) + IF(lwp) WRITE(numout,*)'sbc_fwb : initial freshwater correction flux = ', emp_corr , 'kg/m2/s' + ! + ENDIF + ! + IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN + a_fwb_b = a_fwb ! time swap + agrif_sum = 0._wp ! set work scalar to 0 + CALL Agrif_step_child(glob_sum_volume_agrif) ! integrate over grids + a_fwb = agrif_sum * rho0 / area + IF ( a_fwb_b == 999._wp ) a_fwb_b = a_fwb + emp_corr = (a_fwb - a_fwb_b) / ( rn_Dt * REAL(kn_fsbc, wp) ) + emp_corr + ENDIF + ! + ENDIF + ELSE ! child grid if any + IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN + emp_corr = Agrif_parent(emp_corr) + ENDIF ENDIF ! + IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes on all grids + emp(:,:) = emp(:,:) + emp_corr * tmask(:,:,1) + qns(:,:) = qns(:,:) - emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction + ENDIF + + IF ( Agrif_Root() ) THEN + ! Output restart information (root grid only) + IF( lrst_oce ) THEN + IF(lwp) WRITE(numout,*) + IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt + IF(lwp) WRITE(numout,*) '~~~~' + CALL iom_rstput( kt, nitrst, numrow, 'a_fwb' , a_fwb ) + CALL iom_rstput( kt, nitrst, numrow, 'emp_corr', emp_corr ) + END IF + ! + IF( kt == nitend .AND. lwp ) THEN + WRITE(numout,*) 'sbc_fwb : freshwater-budget at the end of simulation (year now) = ', emp_corr , 'kg/m2/s' + END IF + END IF +#endif + ! outputs + IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN + IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) ) + IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -emp_corr * tmask(:,:,1) ) + ENDIF + ! CASE ( 2 ) !== fw adjustment based on fw budget at the end of the previous year ==! ! simulation is supposed to start 1st of January - IF( kt == nit000 ) THEN ! initialisation - ! ! set the fw adjustment (a_fwb) - IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb_b', ldstop = .FALSE. ) > 0 & ! as read from restart file - & .AND. iom_varid( numror, 'a_fwb', ldstop = .FALSE. ) > 0 ) THEN - IF(lwp) WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file' - CALL iom_get( numror, 'a_fwb_b', a_fwb_b ) - CALL iom_get( numror, 'a_fwb' , a_fwb ) + IF ( Agrif_Root() ) THEN + IF( kt == nit000 ) THEN ! initialisation + ! ! + IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb', ldstop = .FALSE. ) > 0 & ! read from restart file + & .AND. iom_varid( numror, 'a_fwb_b', ldstop = .FALSE. ) > 0 & + & .AND. iom_varid( numror, 'emp_corr', ldstop = .FALSE. ) > 0 ) THEN + IF(lwp) WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file' + CALL iom_get( numror, 'a_fwb' , a_fwb ) + CALL iom_get( numror, 'a_fwb_b' , a_fwb_b ) + CALL iom_get( numror, 'emp_corr' , emp_corr ) + ELSE ! as specified in namelist + IF(lwp) WRITE(numout,*) 'sbc_fwb : setting freshwater-budget from namelist rn_fwb0' + emp_corr = rn_fwb0 + a_fwb = 999._wp + a_fwb_b = 999._wp + END IF ! - a_fwb_ini = a_fwb_b - ELSE ! as specified in namelist - IF(lwp) WRITE(numout,*) 'sbc_fwb : setting freshwater-budget from namelist rn_fwb0' - a_fwb = rn_fwb0 - a_fwb_b = 0._wp ! used only the first year then it is replaced by a_fwb_ini + IF(lwp) WRITE(numout,*) + IF(lwp) WRITE(numout,*)'sbc_fwb : initial freshwater correction flux = ', emp_corr , 'kg/m2/s' ! - a_fwb_ini = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) & - & * rho0 / ( area * rday * REAL(nyear_len(1), wp) ) - END IF - ! - IF(lwp) WRITE(numout,*) - IF(lwp) WRITE(numout,*)'sbc_fwb : freshwater-budget at the end of previous year = ', a_fwb , 'kg/m2/s' - IF(lwp) WRITE(numout,*)' freshwater-budget at initial state = ', a_fwb_ini, 'kg/m2/s' + ENDIF ! - ELSE ! at the end of year n: ikty = nyear_len(1) * 86400 / NINT(rn_Dt) - IF( MOD( kt, ikty ) == 0 ) THEN ! Update a_fwb at the last time step of a year - ! It should be the first time step of a year MOD(kt-1,ikty) but then the restart would be wrong - ! Hence, we make a small error here but the code is restartable - a_fwb_b = a_fwb_ini + IF( MOD( kt-1, ikty ) == 0 ) THEN ! Update a_fwb at the last time step of a year + a_fwb_b = a_fwb ! mean sea level taking into account ice+snow - a_fwb = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) - a_fwb = a_fwb * rho0 / ( area * rday * REAL(nyear_len(1), wp) ) ! convert in kg/m2/s +#if defined key_agrif + agrif_sum = 0.0_wp ! set work scalar to 0 + CALL Agrif_step_child(glob_sum_volume_agrif) ! integrate over grids + a_fwb = agrif_sum +#else + a_fwb = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass_b(:,:) * r1_rho0 ) ) +#endif + a_fwb = a_fwb * rho0 / area + ! + ! Special case if less than a year has been performed: + ! hence namelist rn_fwb0 still rules + IF ( a_fwb_b == 999._wp ) a_fwb_b = a_fwb + ! + emp_corr = ( a_fwb - a_fwb_b ) / ( rday * REAL(nyear_len(0), wp) ) + emp_corr + IF(lwp) WRITE(numout,*) + IF(lwp) WRITE(numout,*)'sbc_fwb : Compute new global mass at step = ', kt + IF(lwp) WRITE(numout,*)'sbc_fwb : New averaged liquid height (ocean + snow + ice) = ', a_fwb * r1_rho0, 'm' + IF(lwp) WRITE(numout,*)'sbc_fwb : Previous averaged liquid height (ocean + snow + ice) = ', a_fwb_b * r1_rho0, 'm' + IF(lwp) WRITE(numout,*)'sbc_fwb : Implied freshwater-budget correction flux = ', emp_corr , 'kg/m2/s' ENDIF - ! +#if defined key_agrif + ELSE ! child grid if any + IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN + emp_corr = Agrif_parent(emp_corr) + ENDIF +#endif ENDIF ! - IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes using previous year budget minus initial state - zcoef = ( a_fwb - a_fwb_b ) - emp(:,:) = emp(:,:) + zcoef * tmask(:,:,1) - qns(:,:) = qns(:,:) - zcoef * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction + IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes + emp(:,:) = emp(:,:) + emp_corr * tmask(:,:,1) + qns(:,:) = qns(:,:) - emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction ! outputs - IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -zcoef * rcp * sst_m(:,:) * tmask(:,:,1) ) - IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -zcoef * tmask(:,:,1) ) - ENDIF - ! Output restart information - IF( lrst_oce ) THEN - IF(lwp) WRITE(numout,*) - IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt - IF(lwp) WRITE(numout,*) '~~~~' - CALL iom_rstput( kt, nitrst, numrow, 'a_fwb_b', a_fwb_b ) - CALL iom_rstput( kt, nitrst, numrow, 'a_fwb', a_fwb ) - END IF - ! - IF( kt == nitend .AND. lwp ) THEN - WRITE(numout,*) 'sbc_fwb : freshwater-budget at the end of simulation (year now) = ', a_fwb , 'kg/m2/s' - WRITE(numout,*) ' freshwater-budget at initial state = ', a_fwb_b, 'kg/m2/s' + IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) ) + IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -emp_corr * tmask(:,:,1) ) ENDIF + + IF ( Agrif_Root() ) THEN + ! Output restart information (root grid only) + IF( lrst_oce ) THEN + IF(lwp) WRITE(numout,*) + IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt + IF(lwp) WRITE(numout,*) '~~~~' + CALL iom_rstput( kt, nitrst, numrow, 'a_fwb', a_fwb ) + CALL iom_rstput( kt, nitrst, numrow, 'a_fwb_b', a_fwb_b ) + CALL iom_rstput( kt, nitrst, numrow, 'emp_corr',emp_corr) + END IF + ! + IF( kt == nitend .AND. lwp ) THEN + IF(lwp) WRITE(numout,*)'sbc_fwb : Previous year averaged liquid height (ocean + snow + ice) = ', a_fwb * r1_rho0, 'm' + IF(lwp) WRITE(numout,*)'sbc_fwb : Previous previous year averaged liquid height (ocean + snow + ice) = ', a_fwb_b * r1_rho0, 'm' + IF(lwp) WRITE(numout,*)'sbc_fwb : freshwater-budget correction flux = ', emp_corr , 'kg/m2/s' + END IF + END IF ! CASE ( 3 ) !== global fwf set to zero and spread out over erp area ==! ! @@ -252,5 +369,30 @@ CONTAINS ! END SUBROUTINE sbc_fwb +#if defined key_agrif + SUBROUTINE glob_sum_area_agrif() + !!--------------------------------------------------------------------- + !! *** compute area with embedded zooms *** + !!---------------------------------------------------------------------- + + agrif_sum = agrif_sum + glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:)) + + END SUBROUTINE glob_sum_area_agrif + + SUBROUTINE glob_sum_volume_agrif() + !!--------------------------------------------------------------------- + !! *** Compute volume with embedded zooms *** + !!---------------------------------------------------------------------- + + IF ( nn_ice==2 ) THEN + agrif_sum = agrif_sum + glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:) * ( ssh(:,:,Kmm_a) + snwice_mass_b(:,:) * r1_rho0 )) + ELSE + agrif_sum = agrif_sum + glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:) * ssh(:,:,Kmm_a) ) + ENDIF + + END SUBROUTINE glob_sum_volume_agrif + +#endif + !!====================================================================== END MODULE sbcfwb diff --git a/src/OCE/SBC/sbcmod.F90 b/src/OCE/SBC/sbcmod.F90 index 29e30de16382c58de2dc09329f82415da1706be3..4b015f7ae984d761235bd1e4baa9ea2cf2a500d0 100644 --- a/src/OCE/SBC/sbcmod.F90 +++ b/src/OCE/SBC/sbcmod.F90 @@ -125,6 +125,10 @@ CONTAINS #endif #if ! defined key_si3 IF( nn_ice == 2 ) nn_ice = 0 ! without key key_si3 you cannot use si3... +#endif +#if defined key_agrif + ! In case of an agrif zoom, the freshwater water budget is determined by parent: + IF (.NOT.Agrif_Root()) nn_fwb = Agrif_Parent(nn_fwb) #endif ! !