MODULE iscplhsb !!====================================================================== !! *** MODULE iscplhsb*** !! Ocean forcing: ice sheet/ocean coupling (conservation) !!===================================================================== !! History : NEMO ! 2015-01 P. Mathiot: original !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! iscpl_alloc : variable allocation !! iscpl_hsb : compute and store the input of heat/salt/volume !! into the system due to the coupling process !! iscpl_div : correction of divergence to keep volume conservation !!---------------------------------------------------------------------- USE dom_oce ! ocean space and time domain USE domwri ! ocean space and time domain USE phycst ! physical constants USE sbc_oce ! surface boundary condition variables 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 domngb ! USE iscplini IMPLICIT NONE PRIVATE PUBLIC iscpl_div PUBLIC iscpl_cons !! * Substitutions # include "domzgr_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_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl) !!---------------------------------------------------------------------- !! *** ROUTINE iscpl_cons *** !! !! ** Purpose : compute input into the system during the coupling step !! compute the correction term !! compute where the correction have to be applied !! !! ** Method : compute tsn*e3t-tsb*e3tb and e3t-e3t_b !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: ptmask_b !! mask before REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pe3t_b !! scale factor before REAL(wp), DIMENSION(:,: ), INTENT(in ) :: psmask_b !! mask before REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx !! corrective flux to have tracer conservation REAL(wp), DIMENSION(:,:,: ), INTENT(out) :: pvol_flx !! corrective flux to have volume conservation REAL(wp), INTENT(in ) :: prdt_iscpl !! coupling period !! INTEGER :: ji, jj, jk !! loop index INTEGER :: jip1, jim1, jjp1, jjm1 !! REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb REAL(wp):: r1_tiscpl REAL(wp):: zjip1_ratio, zjim1_ratio, zjjp1_ratio, zjjm1_ratio !! REAL(wp), DIMENSION(:,: ), POINTER :: zde3t REAL(wp), DIMENSION(:,: ), POINTER :: zssh0 REAL(wp), DIMENSION(:,:,: ), POINTER :: ztmp3d ! REAL(wp), DIMENSION(: ), ALLOCATABLE :: zlon, zlat REAL(wp), DIMENSION(: ), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal INTEGER , DIMENSION(: ), ALLOCATABLE :: ixpts, iypts, izpts, vnpts INTEGER :: jpts, npts CALL wrk_alloc(jpi,jpj,jpk, ztmp3d ) CALL wrk_alloc(jpi,jpj, zde3t ) CALL wrk_alloc(jpi,jpj, zssh0 ) ! get unbalance (volume heat and salt) ! initialisation zde3t (:,:) = 0.0_wp pvol_flx(:,:,: ) = 0.0_wp pts_flx (:,:,:,:) = 0.0_wp zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total volume correction 0 = ',zsum zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total heat correction 0 = ',zsum zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total salt correction 0 = ',zsum ! mask tsn and tsb (should be useless) tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:); tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:); ! diagnose non conservation of heat, salt and volume r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt) zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) IF ( lk_vvl ) zssh0 = 0.0_wp DO jk = 1,jpk-1 DO ji = 2,jpi-1 DO jj = 2,jpj-1 ! volume differences zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk); ! shh changes IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) zssh0(ji,jj) = 0._wp END IF ! ocean cell now ! case where we open, enlarge or thin a cell : pvol_flx(ji,jj,jk) = zde3t(ji,jj) * r1_tiscpl pts_flx (ji,jj,jk,jp_sal)= tsn(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl pts_flx (ji,jj,jk,jp_tem)= tsn(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl END DO END DO END DO ! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0 PRINT *, 'test ', narea, SUM(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total volume correction 1 = ',zsum zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total heat correction 1 = ',zsum zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total salt correction 1 = ',zsum zssh0(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) IF ( lk_vvl ) zssh0 = 0.0_wp DO jk = 1,jpk-1 DO ji = 2,jpi-1 DO jj = 2,jpj-1 ! volume differences zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk); ! shh changes IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) zssh0(ji,jj) = 0._wp END IF ! ocean cell before and mask cell now IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN ! case where we close a cell and adjacent cell open pvol_flx(ji,jj,jk) = zde3t(ji,jj) * r1_tiscpl pts_flx (ji,jj,jk,jp_sal)= tsb(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl pts_flx (ji,jj,jk,jp_tem)= tsb(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ; zsum = e12t(ji ,jjp1) * tmask(ji ,jjp1,jk) + e12t(ji ,jjm1) * tmask(ji ,jjm1,jk) & & + e12t(jim1,jj ) * tmask(jim1,jj ,jk) + e12t(jip1,jj ) * tmask(jip1,jj ,jk) IF ( zsum .NE. 0._wp ) THEN zjip1_ratio = e12t(jip1,jj ) * tmask(jip1,jj ,jk) / zsum zjim1_ratio = e12t(jim1,jj ) * tmask(jim1,jj ,jk) / zsum zjjp1_ratio = e12t(ji ,jjp1) * tmask(ji ,jjp1,jk) / zsum zjjm1_ratio = e12t(ji ,jjm1) * tmask(ji ,jjm1,jk) / zsum pvol_flx(ji ,jjp1,jk ) = pvol_flx(ji ,jjp1,jk ) + pvol_flx(ji,jj,jk ) * zjjp1_ratio pvol_flx(ji ,jjm1,jk ) = pvol_flx(ji ,jjm1,jk ) + pvol_flx(ji,jj,jk ) * zjjm1_ratio pvol_flx(jip1,jj ,jk ) = pvol_flx(jip1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjip1_ratio pvol_flx(jim1,jj ,jk ) = pvol_flx(jim1,jj ,jk ) + pvol_flx(ji,jj,jk ) * zjim1_ratio pts_flx (ji ,jjp1,jk,jp_sal) = pts_flx (ji ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio pts_flx (ji ,jjm1,jk,jp_sal) = pts_flx (ji ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio pts_flx (jip1,jj ,jk,jp_sal) = pts_flx (jip1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio pts_flx (jim1,jj ,jk,jp_sal) = pts_flx (jim1,jj ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio pts_flx (ji ,jjp1,jk,jp_tem) = pts_flx (ji ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio pts_flx (ji ,jjm1,jk,jp_tem) = pts_flx (ji ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio pts_flx (jip1,jj ,jk,jp_tem) = pts_flx (jip1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio pts_flx (jim1,jj ,jk,jp_tem) = pts_flx (jim1,jj ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio ! set to 0 the cell we distributed over neigbourg cells pvol_flx(ji,jj,jk ) = 0._wp pts_flx (ji,jj,jk,jp_sal) = 0._wp pts_flx (ji,jj,jk,jp_tem) = 0._wp ELSE IF (zsum .EQ. 0._wp ) THEN ! case where we close a cell and no adjacent cell open ! check if the cell beneath is wet IF ( tmask(ji,jj,jk+1) .EQ. 1._wp ) THEN pvol_flx(ji,jj,jk+1) = pvol_flx(ji,jj,jk+1) + pvol_flx(ji,jj,jk) pts_flx (ji,jj,jk+1,jp_sal)= pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal) pts_flx (ji,jj,jk+1,jp_tem)= pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem) ! set to 0 the cell we distributed over neigbourg cells pvol_flx(ji,jj,jk ) = 0._wp pts_flx (ji,jj,jk,jp_sal) = 0._wp pts_flx (ji,jj,jk,jp_tem) = 0._wp ELSE ! case no adjacent cell on the horizontal and on the vertical PRINT *, 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal' PRINT *, ' ',mig(ji),' ',mjg(jj),' ',jk PRINT *, ' ',ji,' ',jj,' ',jk,' ',narea PRINT *, ' we are now looking for the closest wet cell on the horizontal ' ! We deal with this points later. END IF END IF END IF END DO END DO END DO zsum = glob_sum_full(pvol_flx(:,:,:) ) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total volume correction 2 = ',zsum zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total heat correction 2 = ',zsum zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt IF (lwp) PRINT *, 'total salt correction 2 = ',zsum ! allocation and initialisation of the list of problematic point ALLOCATE(vnpts(jpnij)) vnpts(:)=0 ! fill narea location with the number of problematic point DO jk = 1,jpk-1 DO ji = 2,jpi-1 DO jj = 2,jpj-1 IF ( ptmask_b(ji,jj,jk) == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & & .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND. tmask(ji,jj,jk+1) == 0 ) THEN vnpts(narea) = vnpts(narea) + 1 END IF END DO END DO END DO ! build array of total problematic point on each cpu (share to each cpu) CALL mpp_max(vnpts,jpnij) ! size of the new variable npts = SUM(vnpts) ! allocation of the coordinates, correction, index vector for the problematic points ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts)) ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20 zcorr_vol(:) = 0.0_wp zcorr_sal(:) = 0.0_wp zcorr_tem(:) = 0.0_wp ! fill new variable jpts = SUM(vnpts(1:narea-1)) DO jk = 1,jpk-1 DO ji = 2,jpi-1 DO jj = 2,jpj-1 IF ( ptmask_b(ji,jj,jk) == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & & .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND. tmask(ji,jj,jk+1) == 0 ) THEN jpts = jpts + 1 ! positioning in the vnpts vector for the area narea PRINT *, 'corrected point ', narea, ji, jj, jk, jpts ixpts(jpts) = ji ; iypts(jpts) = jj ; izpts(jpts) = jk zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj) zcorr_vol(jpts) = pvol_flx(ji,jj,jk) zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal) zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem) ! set flx to 0 (safer) pvol_flx(ji,jj,jk ) = 0.0_wp pts_flx (ji,jj,jk,jp_sal) = 0.0_wp pts_flx (ji,jj,jk,jp_tem) = 0.0_wp PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt END IF END DO END DO END DO ! build array of total problematic point on each cpu (share to each cpu) CALL mpp_max(zlat ,npts) CALL mpp_max(zlon ,npts) CALL mpp_max(izpts,npts) ! put correction term in the closest cell PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts DO jpts = 1,npts CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts)) PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts) DO jj = mj0(iypts(jpts)),mj1(iypts(jpts)) DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts)) jk = izpts(jpts) pvol_flx(ji,jj,jk) = pvol_flx(ji,jj,jk ) + zcorr_vol(jpts) pts_flx (ji,jj,jk,jp_sal) = pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts) pts_flx (ji,jj,jk,jp_tem) = pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts) END DO END DO END DO ! deallocate variables DEALLOCATE(vnpts) DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat) ! add contribution store on the hallo (lbclnk remove one of the contribution) pvol_flx(:,:,: ) = pvol_flx(:,:,: ) * tmask(:,:,:) pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:) pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:) CALL lbc_sum(pvol_flx(:,:,: ),'T',1.) CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.) CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.) ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!! zsumn = glob_sum ( fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt ztmp3d(:,:,:) = 0.0 ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:) zsumb = glob_sum_full(ztmp3d) zsum = glob_sum ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt) IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum ! CHECK salt zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt zsumb = glob_sum( tsb(:,:,:,jp_sal) * pe3t_b(:,:,:) * ptmask_b(:,:,:)) zsum = glob_sum( pts_flx(:,:,:,jp_sal)*rn_fiscpl * rn_rdt) IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum ! CHECK heat zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * tmask (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt zsumb = glob_sum( tsb(:,:,:,jp_tem) * pe3t_b(:,:,:) * ptmask_b(:,:,:)) zsum = glob_sum( pts_flx(:,:,:,jp_tem)*rn_fiscpl * rn_rdt) IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum !! CALL wrk_dealloc(jpi,jpj,jpk, ztmp3d ) CALL wrk_dealloc(jpi,jpj, zde3t ) CALL wrk_dealloc(jpi,jpj, zssh0 ) END SUBROUTINE iscpl_cons SUBROUTINE iscpl_div( phdivn ) !!---------------------------------------------------------------------- !! *** ROUTINE iscpl_div *** !! !! ** Purpose : update the horizontal divergenc !! !! ** 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 !! INTEGER :: ji, jj, jk ! dummy loop indices !!---------------------------------------------------------------------- ! DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk) END DO END DO END DO ! END SUBROUTINE iscpl_div END MODULE iscplhsb