#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES */ #undef DECAL_FEEDBACK_2D /* SEPARATION of INTERFACES (Barotropic mode) */ #undef VOL_REFLUX /* VOLUME REFLUXING*/ MODULE agrif_oce_update !!====================================================================== !! *** MODULE agrif_oce_interp *** !! AGRIF: update package for the ocean dynamics (OPA) !!====================================================================== !! History : 2.0 ! 2002-06 (L. Debreu) Original code !! 3.2 ! 2009-04 (R. Benshila) !! 3.6 ! 2014-09 (R. Benshila) !!---------------------------------------------------------------------- #if defined key_agrif !!---------------------------------------------------------------------- !! 'key_agrif' AGRIF zoom !!---------------------------------------------------------------------- USE par_oce USE oce USE dom_oce USE zdf_oce ! vertical physics: ocean variables USE agrif_oce ! USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE domvvl ! Need interpolation routines USE vremap ! Vertical remapping USE lbclnk IMPLICIT NONE PRIVATE PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh PUBLIC Update_Scales !!---------------------------------------------------------------------- !! NEMO/NST 4.0 , NEMO Consortium (2018) !! $Id$ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE Agrif_Update_Tra( ) !!---------------------------------------------------------------------- !! *** ROUTINE Agrif_Update_Tra *** !!---------------------------------------------------------------------- ! IF (Agrif_Root()) RETURN ! IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers from grid Number',Agrif_Fixed() #if defined key_vertical ! Effect of this has to be carrefully checked ! depending on what the nesting tools ensure for ! volume conservation: Agrif_UseSpecialValueInUpdate = .FALSE. #else Agrif_UseSpecialValueInUpdate = .TRUE. #endif Agrif_SpecialValueFineGrid = 0._wp ! # if ! defined DECAL_FEEDBACK CALL Agrif_Update_Variable(tsn_id, procname=updateTS) ! near boundary update: ! CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS) # else CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS) ! near boundary update: ! CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS) # endif ! Agrif_UseSpecialValueInUpdate = .FALSE. ! ! END SUBROUTINE Agrif_Update_Tra SUBROUTINE Agrif_Update_Dyn( ) !!---------------------------------------------------------------------- !! *** ROUTINE Agrif_Update_Dyn *** !!---------------------------------------------------------------------- ! IF (Agrif_Root()) RETURN ! IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed() Agrif_UseSpecialValueInUpdate = .FALSE. Agrif_SpecialValueFineGrid = 0. use_sign_north = .TRUE. sign_north = -1. ! # if ! defined DECAL_FEEDBACK CALL Agrif_Update_Variable(un_update_id,procname = updateU) CALL Agrif_Update_Variable(vn_update_id,procname = updateV) ! near boundary update: ! CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU) ! CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV) # else CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU) CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV) ! near boundary update: ! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU) ! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) # endif # if ! defined DECAL_FEEDBACK_2D CALL Agrif_Update_Variable(e1u_id,procname = updateU2d) CALL Agrif_Update_Variable(e2v_id,procname = updateV2d) # else CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d) CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d) # endif ! # if ! defined DECAL_FEEDBACK_2D ! Account for updated thicknesses at boundary edges IF (.NOT.ln_linssh) THEN ! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy) ! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy) ENDIF # endif ! IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN ! Update time integrated transports # if ! defined DECAL_FEEDBACK_2D CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) # else CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b) CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b) # endif END IF ! use_sign_north = .FALSE. ! END SUBROUTINE Agrif_Update_Dyn SUBROUTINE Agrif_Update_ssh( ) !!--------------------------------------------- !! *** ROUTINE Agrif_Update_ssh *** !!--------------------------------------------- ! IF (Agrif_Root()) RETURN ! Agrif_UseSpecialValueInUpdate = .TRUE. Agrif_SpecialValueFineGrid = 0. # if ! defined DECAL_FEEDBACK_2D CALL Agrif_Update_Variable(sshn_id,procname = updateSSH) # else CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH) # endif ! Agrif_UseSpecialValueInUpdate = .FALSE. ! # if defined VOL_REFLUX IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN use_sign_north = .TRUE. sign_north = -1. ! Refluxing on ssh: # if defined DECAL_FEEDBACK_2D CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) # else CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu) CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv) # endif use_sign_north = .FALSE. END IF # endif ! END SUBROUTINE Agrif_Update_ssh SUBROUTINE Agrif_Update_Tke( ) !!--------------------------------------------- !! *** ROUTINE Agrif_Update_Tke *** !!--------------------------------------------- !! ! IF (Agrif_Root()) RETURN ! Agrif_UseSpecialValueInUpdate = .TRUE. Agrif_SpecialValueFineGrid = 0. CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN ) CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT ) CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM ) Agrif_UseSpecialValueInUpdate = .FALSE. END SUBROUTINE Agrif_Update_Tke SUBROUTINE Agrif_Update_vvl( ) !!--------------------------------------------- !! *** ROUTINE Agrif_Update_vvl *** !!--------------------------------------------- ! IF (Agrif_Root()) RETURN ! IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() ! Agrif_UseSpecialValueInUpdate = .TRUE. Agrif_SpecialValueFineGrid = 0. ! ! No interface separation here, update vertical grid at T points ! everywhere over the overlapping regions (one account for refluxing in that case): CALL Agrif_Update_Variable(e3t_id, procname=updatee3t) ! Agrif_UseSpecialValueInUpdate = .FALSE. ! CALL Agrif_ChildGrid_To_ParentGrid() CALL dom_vvl_update_UVF CALL Agrif_ParentGrid_To_ChildGrid() ! END SUBROUTINE Agrif_Update_vvl SUBROUTINE dom_vvl_update_UVF !!--------------------------------------------- !! *** ROUTINE dom_vvl_update_UVF *** !!--------------------------------------------- !! INTEGER :: jk REAL(wp):: zcoef !!--------------------------------------------- IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', & & Agrif_Fixed(), 'Step', Agrif_Nb_Step() ! Save "old" scale factor (prior update) for subsequent asselin correction ! of prognostic variables ! ----------------------- ! e3u(:,:,:,Krhs_a) = e3u(:,:,:,Kmm_a) e3v(:,:,:,Krhs_a) = e3v(:,:,:,Kmm_a) ! uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a) ! vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a) hu(:,:,Krhs_a) = hu(:,:,Kmm_a) hv(:,:,Krhs_a) = hv(:,:,Kmm_a) ! 1) NOW fields !-------------- ! Vertical scale factor interpolations ! ------------------------------------ CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3u(:,:,:,Kmm_a) , 'U' ) CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3v(:,:,:,Kmm_a) , 'V' ) CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3f(:,:,:) , 'F' ) CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3uw(:,:,:,Kmm_a), 'UW' ) CALL dom_vvl_interpol( e3v(:,:,:,Kmm_a), e3vw(:,:,:,Kmm_a), 'VW' ) ! Update total depths: ! -------------------- hu(:,:,Kmm_a) = 0._wp ! Ocean depth at U-points hv(:,:,Kmm_a) = 0._wp ! Ocean depth at V-points DO jk = 1, jpkm1 hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk) hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk) END DO ! ! Inverse of the local depth r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) ) r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) ) ! 2) BEFORE fields: !------------------ IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN ! ! Vertical scale factor interpolations ! ------------------------------------ CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3u(:,:,:,Kbb_a), 'U' ) CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3v(:,:,:,Kbb_a), 'V' ) CALL dom_vvl_interpol( e3u(:,:,:,Kbb_a), e3uw(:,:,:,Kbb_a), 'UW' ) CALL dom_vvl_interpol( e3v(:,:,:,Kbb_a), e3vw(:,:,:,Kbb_a), 'VW' ) ! Update total depths: ! -------------------- hu(:,:,Kbb_a) = 0._wp ! Ocean depth at U-points hv(:,:,Kbb_a) = 0._wp ! Ocean depth at V-points DO jk = 1, jpkm1 hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk) hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk) END DO ! ! Inverse of the local depth r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) ) r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) ) ENDIF ! END SUBROUTINE dom_vvl_update_UVF #if defined key_vertical SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE updateT *** !!--------------------------------------------- INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL, INTENT(in) :: before !! INTEGER :: ji,jj,jk,jn INTEGER :: N_in, N_out REAL(wp) :: ztb, ztnu, ztno REAL(wp) :: h_in(k1:k2) REAL(wp) :: h_out(1:jpk) REAL(wp) :: tabin(k1:k2,1:jpts) REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: tabres_child !!--------------------------------------------- ! IF (before) THEN !jc_alt ! AGRIF_SpecialValue = -999._wp DO jn = n1,n2-1 DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 !jc_alt ! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & ! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) END DO END DO END DO END DO DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 !jc_alt ! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & ! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) END DO END DO END DO ELSE tabres_child(:,:,:,:) = 0._wp AGRIF_SpecialValue = 0._wp DO jj=j1,j2 DO ji=i1,i2 N_in = 0 DO jk=k1,k2 !k2 = jpk of child grid ! jc_alt ! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT N_in = N_in + 1 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) h_in(N_in) = tabres(ji,jj,jk,n2) ENDDO N_out = 0 DO jk=1,jpk ! jpk of parent grid IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF N_out = N_out + 1 h_out(N_out) = e3t(ji,jj,jk,Kmm_a) ENDDO IF (N_in*N_out > 0) THEN !Remove this? CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) ENDIF ENDDO ENDDO IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part DO jn = 1,jpts DO jk = 1, jpkm1 DO jj = j1, j2 DO ji = i1, i2 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) & & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) ENDIF END DO END DO END DO END DO ENDIF DO jn = 1,jpts DO jk = 1, jpkm1 DO jj = j1, j2 DO ji = i1, i2 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) END IF END DO END DO END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a) ENDIF ENDIF ! END SUBROUTINE updateTS # else SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) !!--------------------------------------------- !! *** ROUTINE updateT *** !!--------------------------------------------- INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL, INTENT(in) :: before !! INTEGER :: ji,jj,jk,jn REAL(wp) :: ztb, ztnu, ztno !!--------------------------------------------- ! IF (before) THEN DO jn = 1,jpts DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 !> jc tmp tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) ! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) !< jc tmp END DO END DO END DO END DO ELSE !> jc tmp DO jn = 1,jpts tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & & * tmask(i1:i2,j1:j2,k1:k2) ENDDO !< jc tmp IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part DO jn = 1,jpts DO jk = k1, k2 DO jj = j1, j2 DO ji = i1, i2 IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used ztnu = tabres(ji,jj,jk,jn) ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) & & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) ENDIF END DO END DO END DO END DO ENDIF DO jn = 1,jpts DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) END IF END DO END DO END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a) ENDIF ! ENDIF ! END SUBROUTINE updateTS # endif # if defined key_vertical SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) !!--------------------------------------------- !! *** ROUTINE updateu *** !!--------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before ! INTEGER :: ji, jj, jk REAL(wp):: zrhoy, zub, zunu, zuno ! VERTICAL REFINEMENT BEGIN REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child REAL(wp) :: h_in(k1:k2) REAL(wp) :: h_out(1:jpk) INTEGER :: N_in, N_out REAL(wp) :: h_diff, excess, thick REAL(wp) :: tabin(k1:k2) ! VERTICAL REFINEMENT END !!--------------------------------------------- ! IF( before ) THEN zrhoy = Agrif_Rhoy() !jc_alt ! AGRIF_SpecialValue = -999._wp DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 !jc_alt ! tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) & ! & + (umask(ji,jj,jk)-1._wp)*999._wp tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) !jc_alt ! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & ! & + (umask(ji,jj,jk)-1._wp)*999._wp tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) END DO END DO END DO ELSE tabres_child(:,:,:) = 0. AGRIF_SpecialValue = 0._wp DO jj=j1,j2 DO ji=i1,i2 N_in = 0 h_in(:) = 0._wp tabin(:) = 0._wp DO jk=k1,k2 !k2=jpk of child grid !jc_alt ! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT IF( tabres(ji,jj,jk,2) == 0.) EXIT N_in = N_in + 1 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) ENDDO N_out = 0 DO jk=1,jpk IF (umask(ji,jj,jk) == 0) EXIT N_out = N_out + 1 h_out(N_out) = e3u(ji,jj,jk,Kmm_a) ENDDO IF (N_in * N_out > 0) THEN h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) excess = 0._wp IF (h_diff < -1.e-4) THEN !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. DO jk=N_in,1,-1 thick = MIN(-1*h_diff, h_in(jk)) excess = excess + tabin(jk)*thick*e2u(ji,jj) tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) h_diff = h_diff + thick IF ( h_diff == 0) THEN N_in = jk h_in(jk) = h_in(jk) - thick EXIT ENDIF ENDDO ENDIF CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) ENDIF ENDDO ENDDO ! DO jk=1,jpk DO jj=j1,j2 DO ji=i1,i2 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a) uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) & & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) ENDIF ! uu(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) END DO END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) ENDIF ! ENDIF ! END SUBROUTINE updateu #else SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) !!--------------------------------------------- !! *** ROUTINE updateu *** !!--------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before ! INTEGER :: ji, jj, jk REAL(wp) :: zrhoy, zub, zunu, zuno !!--------------------------------------------- ! IF( before ) THEN zrhoy = Agrif_Rhoy() DO jk = k1, k2 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a) END DO ELSE DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) ! IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) zunu = tabres(ji,jj,jk,1) uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) & & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) ENDIF ! uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a) END DO END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN uu(i1:i2,j1:j2,k1:k2,Kbb_a) = uu(i1:i2,j1:j2,k1:k2,Kmm_a) ENDIF ! ENDIF ! END SUBROUTINE updateu # endif SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) !!--------------------------------------------- !! *** ROUTINE correct_u_bdy *** !!--------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before INTEGER , INTENT(in) :: nb, ndir !! LOGICAL :: western_side, eastern_side ! INTEGER :: jj, jk REAL(wp) :: zcor !!--------------------------------------------- ! IF( .NOT.before ) THEN ! western_side = (nb == 1).AND.(ndir == 1) eastern_side = (nb == 1).AND.(ndir == 2) ! IF (western_side) THEN DO jj=j1,j2 zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a) uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor DO jk=1,jpkm1 uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk) END DO END DO ENDIF ! IF (eastern_side) THEN DO jj=j1,j2 zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a) uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor DO jk=1,jpkm1 uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk) END DO END DO ENDIF ! ENDIF ! END SUBROUTINE correct_u_bdy # if defined key_vertical SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) !!--------------------------------------------- !! *** ROUTINE updatev *** !!--------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before ! INTEGER :: ji, jj, jk REAL(wp) :: zrhox, zvb, zvnu, zvno ! VERTICAL REFINEMENT BEGIN REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child REAL(wp) :: h_in(k1:k2) REAL(wp) :: h_out(1:jpk) INTEGER :: N_in, N_out REAL(wp) :: h_diff, excess, thick REAL(wp) :: tabin(k1:k2) ! VERTICAL REFINEMENT END !!--------------------------------------------- ! IF( before ) THEN zrhox = Agrif_Rhox() !jc_alt ! AGRIF_SpecialValue = -999._wp DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 !jc_alt ! tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) & ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) !jc_alt ! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) & ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) END DO END DO END DO ELSE tabres_child(:,:,:) = 0. AGRIF_SpecialValue = 0._wp DO jj=j1,j2 DO ji=i1,i2 N_in = 0 DO jk=k1,k2 !jc_alt ! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT IF (tabres(ji,jj,jk,2) == 0) EXIT N_in = N_in + 1 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) ENDDO N_out = 0 DO jk=1,jpk IF (vmask(ji,jj,jk) == 0) EXIT N_out = N_out + 1 h_out(N_out) = e3v(ji,jj,jk,Kmm_a) ENDDO IF (N_in * N_out > 0) THEN h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) excess = 0._wp IF (h_diff < -1.e-4) then !Even if bathy at T points match it's possible for the V points to be deeper in the child grid. !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. DO jk=N_in,1,-1 thick = MIN(-1*h_diff, h_in(jk)) excess = excess + tabin(jk)*thick*e2u(ji,jj) tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) h_diff = h_diff + thick IF ( h_diff == 0) THEN N_in = jk h_in(jk) = h_in(jk) - thick EXIT ENDIF ENDDO ENDIF CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) ENDIF ENDDO ENDDO ! DO jk=1,jpkm1 DO jj=j1,j2 DO ji=i1,i2 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a) vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) & & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) ENDIF ! vv(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) END DO END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) ENDIF ! ENDIF ! END SUBROUTINE updatev # else SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) !!--------------------------------------------- !! *** ROUTINE updatev *** !!--------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before ! INTEGER :: ji, jj, jk REAL(wp) :: zrhox, zvb, zvnu, zvno !!--------------------------------------------- ! IF (before) THEN zrhox = Agrif_Rhox() DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) END DO END DO END DO ELSE DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) ! IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) zvnu = tabres(ji,jj,jk,1) vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) & & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) ENDIF ! vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a) END DO END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN vv(i1:i2,j1:j2,k1:k2,Kbb_a) = vv(i1:i2,j1:j2,k1:k2,Kmm_a) ENDIF ! ENDIF ! END SUBROUTINE updatev # endif SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) !!--------------------------------------------- !! *** ROUTINE correct_v_bdy *** !!--------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before INTEGER , INTENT(in) :: nb, ndir !! LOGICAL :: southern_side, northern_side ! INTEGER :: ji, jk REAL(wp) :: zcor !!--------------------------------------------- ! IF( .NOT.before ) THEN ! southern_side = (nb == 2).AND.(ndir == 1) northern_side = (nb == 2).AND.(ndir == 2) ! IF (southern_side) THEN DO ji=i1,i2 zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a) vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor DO jk=1,jpkm1 vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk) END DO END DO ENDIF ! IF (northern_side) THEN DO ji=i1,i2 zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a) vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor DO jk=1,jpkm1 vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk) END DO END DO ENDIF ! ENDIF ! END SUBROUTINE correct_v_bdy SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE updateu2d *** !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before !! INTEGER :: ji, jj, jk REAL(wp) :: zrhoy REAL(wp) :: zcorr !!--------------------------------------------- ! IF( before ) THEN zrhoy = Agrif_Rhoy() DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj) END DO END DO ELSE DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj) ! ! Update "now" 3d velocities: spgu(ji,jj) = 0._wp DO jk=1,jpkm1 spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) END DO ! zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) DO jk=1,jpkm1 uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk) END DO ! ! Update barotropic velocities: IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a) uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + rn_atfp * zcorr * umask(ji,jj,1) END IF ENDIF uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) ! ! Correct "before" velocities to hold correct bt component: spgu(ji,jj) = 0.e0 DO jk=1,jpkm1 spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) END DO ! zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a) DO jk=1,jpkm1 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk) END DO ! END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN uu_b(i1:i2,j1:j2,Kbb_a) = uu_b(i1:i2,j1:j2,Kmm_a) ENDIF ENDIF ! END SUBROUTINE updateu2d SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE updatev2d *** !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before ! INTEGER :: ji, jj, jk REAL(wp) :: zrhox, zcorr !!---------------------------------------------------------------------- ! IF( before ) THEN zrhox = Agrif_Rhox() DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj) END DO END DO ELSE DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj) ! ! Update "now" 3d velocities: spgv(ji,jj) = 0.e0 DO jk=1,jpkm1 spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) END DO ! zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) DO jk=1,jpkm1 vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk) END DO ! ! Update barotropic velocities: IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a) vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + rn_atfp * zcorr * vmask(ji,jj,1) END IF ENDIF vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) ! ! Correct "before" velocities to hold correct bt component: spgv(ji,jj) = 0.e0 DO jk=1,jpkm1 spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) END DO ! zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a) DO jk=1,jpkm1 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk) END DO ! END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN vv_b(i1:i2,j1:j2,Kbb_a) = vv_b(i1:i2,j1:j2,Kmm_a) ENDIF ! ENDIF ! END SUBROUTINE updatev2d SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE updateSSH *** !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before !! INTEGER :: ji, jj !!---------------------------------------------------------------------- ! IF( before ) THEN DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = ssh(ji,jj,Kmm_a) END DO END DO ELSE IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN DO jj=j1,j2 DO ji=i1,i2 ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) & & + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1) END DO END DO ENDIF ! DO jj=j1,j2 DO ji=i1,i2 ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1) END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN ssh(i1:i2,j1:j2,Kbb_a) = ssh(i1:i2,j1:j2,Kmm_a) ENDIF ! ENDIF ! END SUBROUTINE updateSSH SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE updateub2b *** !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: i1, i2, j1, j2 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres LOGICAL , INTENT(in) :: before !! INTEGER :: ji, jj REAL(wp) :: zrhoy, za1, zcor !!--------------------------------------------- ! IF (before) THEN zrhoy = Agrif_Rhoy() DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj) END DO END DO tabres = zrhoy * tabres ELSE ! tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2) ! za1 = 1._wp / REAL(Agrif_rhot(), wp) DO jj=j1,j2 DO ji=i1,i2 zcor=tabres(ji,jj) - ub2_b(ji,jj) ! Update time integrated fluxes also in case of multiply nested grids: ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor ! Update corrective fluxes: un_bf(ji,jj) = un_bf(ji,jj) + zcor ! Update half step back fluxes: ub2_b(ji,jj) = tabres(ji,jj) END DO END DO ENDIF ! END SUBROUTINE updateub2b SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir ) !!--------------------------------------------- !! *** ROUTINE reflux_sshu *** !!--------------------------------------------- INTEGER, INTENT(in) :: i1, i2, j1, j2 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres LOGICAL, INTENT(in) :: before INTEGER, INTENT(in) :: nb, ndir !! LOGICAL :: western_side, eastern_side INTEGER :: ji, jj REAL(wp) :: zrhoy, za1, zcor !!--------------------------------------------- ! IF (before) THEN zrhoy = Agrif_Rhoy() DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj) END DO END DO tabres = zrhoy * tabres ELSE ! tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2) ! western_side = (nb == 1).AND.(ndir == 1) eastern_side = (nb == 1).AND.(ndir == 2) ! IF (western_side) THEN DO jj=j1,j2 zcor = rn_Dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + rn_atfp * zcor END DO ENDIF IF (eastern_side) THEN DO jj=j1,j2 zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + rn_atfp * zcor END DO ENDIF ! ENDIF ! END SUBROUTINE reflux_sshu SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE updatevb2b *** !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before !! INTEGER :: ji, jj REAL(wp) :: zrhox, za1, zcor !!--------------------------------------------- ! IF( before ) THEN zrhox = Agrif_Rhox() DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) END DO END DO tabres = zrhox * tabres ELSE ! tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2) ! za1 = 1._wp / REAL(Agrif_rhot(), wp) DO jj=j1,j2 DO ji=i1,i2 zcor=tabres(ji,jj) - vb2_b(ji,jj) ! Update time integrated fluxes also in case of multiply nested grids: vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor ! Update corrective fluxes: vn_bf(ji,jj) = vn_bf(ji,jj) + zcor ! Update half step back fluxes: vb2_b(ji,jj) = tabres(ji,jj) END DO END DO ENDIF ! END SUBROUTINE updatevb2b SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir ) !!--------------------------------------------- !! *** ROUTINE reflux_sshv *** !!--------------------------------------------- INTEGER, INTENT(in) :: i1, i2, j1, j2 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres LOGICAL, INTENT(in) :: before INTEGER, INTENT(in) :: nb, ndir !! LOGICAL :: southern_side, northern_side INTEGER :: ji, jj REAL(wp) :: zrhox, za1, zcor !!--------------------------------------------- ! IF (before) THEN zrhox = Agrif_Rhox() DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) END DO END DO tabres = zrhox * tabres ELSE ! tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2) ! southern_side = (nb == 2).AND.(ndir == 1) northern_side = (nb == 2).AND.(ndir == 2) ! IF (southern_side) THEN DO ji=i1,i2 zcor = rn_Dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1)) ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + rn_atfp * zcor END DO ENDIF IF (northern_side) THEN DO ji=i1,i2 zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2)) ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + rn_atfp * zcor END DO ENDIF ! ENDIF ! END SUBROUTINE reflux_sshv SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before ) ! ! ====>>>>>>>>>> currently not used ! !!---------------------------------------------------------------------- !! *** ROUTINE updateT *** !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres LOGICAL , INTENT(in ) :: before !! INTEGER :: ji,jj,jk REAL(wp) :: ztemp !!---------------------------------------------------------------------- IF (before) THEN DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk) tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk) END DO END DO END DO tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy() tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox() tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy() ELSE DO jk=k1,k2 DO jj=j1,j2 DO ji=i1,i2 IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk) print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk) ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3))) print *,'CORR = ',ztemp-1. print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, & tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp END IF END DO END DO END DO ENDIF ! END SUBROUTINE update_scales SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE updateen *** !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab LOGICAL , INTENT(in ) :: before !!---------------------------------------------------------------------- ! IF( before ) THEN ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2) ELSE en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) ENDIF ! END SUBROUTINE updateEN SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before ) !!---------------------------------------------------------------------- !! *** ROUTINE updateavt *** !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab LOGICAL , INTENT(in ) :: before !!---------------------------------------------------------------------- ! IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2) ELSE ; avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) ENDIF ! END SUBROUTINE updateAVT SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before ) !!--------------------------------------------- !! *** ROUTINE updateavm *** !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab LOGICAL , INTENT(in ) :: before !!---------------------------------------------------------------------- ! IF( before ) THEN ; ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2) ELSE ; avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) ENDIF ! END SUBROUTINE updateAVM SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before ) !!--------------------------------------------- !! *** ROUTINE updatee3t *** !!--------------------------------------------- REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2 LOGICAL, INTENT(in) :: before ! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab INTEGER :: ji,jj,jk REAL(wp) :: zcoef !!--------------------------------------------- ! IF (.NOT.before) THEN ! ALLOCATE(ptab(i1:i2,j1:j2,1:jpk)) ! ! Update e3t from ssh (z* case only) DO jk = 1, jpkm1 DO jj=j1,j2 DO ji=i1,i2 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) & & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj))) END DO END DO END DO ! ! 1) Updates at BEFORE time step: ! ------------------------------- ! ! Save "old" scale factor (prior update) for subsequent asselin correction ! of prognostic variables e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) ! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace... ! hdiv(i1:i2,j1:j2,1:jpkm1) = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a) IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN DO jk = 1, jpkm1 DO jj=j1,j2 DO ji=i1,i2 e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) & & + rn_atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) ) END DO END DO END DO ! e3w (i1:i2,j1:j2,1,Kbb_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kbb_a) - e3t_0(i1:i2,j1:j2,1) gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a) ! DO jk = 2, jpk DO jj = j1,j2 DO ji = i1,i2 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * & & ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) & & + 0.5_wp * tmask(ji,jj,jk) * & & ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) ) gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a) gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) & & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a)) END DO END DO END DO ! ENDIF ! ! 2) Updates at NOW time step: ! ---------------------------- ! ! Update vertical scale factor at T-points: e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1) ! ! Update total depth: ht(i1:i2,j1:j2) = 0._wp DO jk = 1, jpkm1 ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk) END DO ! ! Update vertical scale factor at W-points and depths: e3w (i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1) gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a) gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm_a) - (ht(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh ! DO jk = 2, jpk DO jj = j1,j2 DO ji = i1,i2 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) e3w(ji,jj,jk,Kmm_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) ) & & + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) ) gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a) gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) & & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a)) gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh END DO END DO END DO ! IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a) e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a) gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a) gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a) ENDIF ! DEALLOCATE(ptab) ENDIF ! END SUBROUTINE updatee3t #else !!---------------------------------------------------------------------- !! Empty module no AGRIF zoom !!---------------------------------------------------------------------- CONTAINS SUBROUTINE agrif_oce_update_empty WRITE(*,*) 'agrif_oce_update : You should not have seen this print! error?' END SUBROUTINE agrif_oce_update_empty #endif !!====================================================================== END MODULE agrif_oce_update