MODULE iscpldiv !!====================================================================== !! *** MODULE iscpldiv*** !! Divergence correction : compute and save the divergence correction due to ice sheet/ocean coupling !!===================================================================== !! History : NEMO ! 2018-06 A. Siahaan: original !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! iscpl_div_corr : compute the divergence correction !! iscpl_div : update hdivn with the divergence correction !! iscpl_div_tra : update tsa with contribution from the divergence correction !!---------------------------------------------------------------------- USE dom_oce ! ocean space and time domain USE oce ! global tra/dyn variable USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE lib_fortran ! MPP library USE wrk_nemo ! Memory allocation USE lbclnk ! USE iscplini IMPLICIT NONE PRIVATE PUBLIC iscpl_div PUBLIC iscpl_div_tra PUBLIC iscpl_div_corr !! !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE iscpl_div( phdivn ) !!---------------------------------------------------------------------- !! *** ROUTINE iscpl_div *** !! !! ** Purpose : update the horizontal divergence !! !! ** Method : !! CAUTION : iscpl is positive (inflow) and expressed in !m/s !! !! ** Action : phdivn increase by the iscpl correction term !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence !! !!---------------------------------------------------------------------- ! IF(lwp) WRITE(numout,*) 'adding the divergence correction contribution at nit000' phdivn(:,:,:) = phdivn(:,:,:) + rhdivdiff(:,:,:) END SUBROUTINE iscpl_div SUBROUTINE iscpl_div_tra( ptsa ) !!---------------------------------------------------------------------- !! *** ROUTINE iscpl_div_tra *** !! !! ** Purpose : update the active tracer tendency !! !! ** Method : !! CAUTION : iscpl is positive (inflow) and expressed in !m/s !! !! ** Action : ptsa increase by the iscpl correction term !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: ptsa ! active tracer tendency !! !!---------------------------------------------------------------------- ! IF (lwp) WRITE(numout,*) 'add tendency tsa due to re-meshing at kt = nit000' ptsa(:,:,:,1) = ptsa(:,:,:,1) + rhdivdiff_trac(:,:,:,1)*tmask(:,:,:) ptsa(:,:,:,2) = ptsa(:,:,:,2) + rhdivdiff_trac(:,:,:,2)*tmask(:,:,:) END SUBROUTINE iscpl_div_tra SUBROUTINE iscpl_div_corr(ptmask_b, pe3t_b, pe3u_b, pe3v_b) !!---------------------------------------------------------------------- !! *** ROUTINE iscpl_div_corr *** !! !! ** Purpose : compute the divergence correction !! !! ** Method : !! compute the difference of hdivn*e3t_n between before and after remeshing !! transfer the missing hdvn*e3t_n from top of isf to the new mikt !! transfer the missing hdvn*e3t_n from bottom of isf to the new mbkt !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: ptmask_b !! mask before REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pe3t_b , pe3u_b , pe3v_b !! scale factor before !! INTEGER :: ji, jj, jk, jtop !! loop index !! REAL(wp):: zvol, zuflux_sum, zvflux_sum REAL(wp), DIMENSION(:,:,: ), POINTER :: zhdiv, zhdiv2 !!---------------------------------------------------------------------- !! allocate variables CALL wrk_alloc(jpi,jpj,jpk, zhdiv, zhdiv2) !-----initialise the divergence and tracer correction------ zhdiv(:,:,:) = 0._wp rhdivdiff(:,:,:) = 0._wp rhdivdiff_trac(:,:,:,:) = 0._wp zhdiv2(:,:,:) = 0._wp !calculate the hdivn*e3tn with the last values before re-meshing (should only be different from the real hdivn for cavity only) !hdivb and rot should remain the same as the previous timestep bcause it will be used in dyn_ldf, so we did not use div_cur to calculate rhdivdiff as we should only use div_cur once for div_cur(nit000) in order not to update hdivb and rot !----------------------------------------------hdivn before remeshing------------------------------------- DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpk ! Horizontal slab zhdiv2(ji,jj,jk) = & ( e2u(ji,jj)*pe3u_b(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj )*pe3u_b(ji-1,jj ,jk) * un(ji-1,jj ,jk) & + e1v(ji,jj)*pe3v_b(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji ,jj-1)*pe3v_b(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) & / ( e1t(ji,jj) * e2t(ji,jj) )*ptmask_b(ji,jj,jk) END DO END DO ENDDO ! velocity with the new mask un(:,:,:)=un(:,:,:)*umask(:,:,:) vn(:,:,:)=vn(:,:,:)*vmask(:,:,:) !calling divcur to update hdivn (ingat: hdivb = hdivn(before re-meshing), yet that hdivn was not using the recent values, as it was not updated at the end of timestep ) !CALL div_cur(nit000) !----------------------------------------------hdivn after remeshing------------------------------------- DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpk ! Horizontal slab zhdiv(ji,jj,jk) = & ( e2u(ji,jj)*fse3u_n(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj )*fse3u_n(ji-1,jj ,jk) * un(ji-1,jj ,jk) & + e1v(ji,jj)*fse3v_n(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji ,jj-1)*fse3v_n(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) & / ( e1t(ji,jj) * e2t(ji,jj) )*tmask(ji,jj,jk) END DO END DO ENDDO zhdiv2(:,:,:) = zhdiv2(:,:,:) - zhdiv(:,:,:) ! non-cavity will be = 0 rhdivdiff(:,:,:) = zhdiv2(:,:,:) !----transfers the divergence of closed cells to the new top cell------ !until this point, rhdivdiff has the same unit as hdivn*e3t DO jj = 2, jpjm1 DO ji = 2, jpim1 jk = mikt(ji,jj) IF (jk > 1 .AND. ssmask(ji,jj) == 1 ) THEN rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + SUM(rhdivdiff(ji,jj,1:jk-1)*ptmask_b(ji,jj,1:jk-1)) ENDIF jk = mbkt(ji,jj) IF (ssmask(ji,jj) == 1 ) THEN rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + SUM(rhdivdiff(ji,jj,jk+1:jpk)*ptmask_b(ji,jj,jk+1:jpk)) ENDIF END DO ENDDO !----cancel all the divergence of closed cells with tmask, then convert it to hdivn by dividing with e3t_n (it is now safe to do this) ! also cancelling totally closed cavity column rhdivdiff(:,:,:) = rhdivdiff(:,:,:)*tmask(:,:,:)/fse3t_n(:,:,:) !tracer correction trend associated with the divergence correction; opposite sign of divergence rhdivdiff_trac(:,:,:,1) = - rhdivdiff(:,:,:) * tsb(:,:,:,1)*tmask(:,:,:) rhdivdiff_trac(:,:,:,2) = - rhdivdiff(:,:,:) * tsb(:,:,:,2)*tmask(:,:,:) CALL lbc_lnk(rhdivdiff, 'T',1._wp) CALL lbc_lnk(rhdivdiff_trac(:,:,:,1), 'T',1._wp) CALL lbc_lnk(rhdivdiff_trac(:,:,:,2), 'T',1._wp) CALL wrk_dealloc(jpi,jpj,jpk, zhdiv, zhdiv2 ) END SUBROUTINE iscpl_div_corr END MODULE iscpldiv