MODULE domvvl !!====================================================================== !! *** MODULE domvvl *** !! Ocean : !!====================================================================== !! History : 9.0 ! 06-06 (B. Levier, L. Marie) original code !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! 'key_vvl' variable volume !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dom_vvl : defined scale factors & depths at each time step !! dom_vvl_ini : defined coefficients to distribute ssh on each layers !! dom_vvl_ssh : defined the ocean sea level at each time step !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE in_out_manager ! I/O manager USE lib_mpp ! distributed memory computing library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE dynspg_oce ! surface pressure gradient variables USE ocesbc ! ocean surface boundary condition USE phycst ! physical constants IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC dom_vvl_ini ! called by dom_init.F90 PUBLIC dom_vvl_ssh ! called by trazdf.F90 PUBLIC dom_vvl ! called by istate.F90 and step.F90 PUBLIC dom_vvl_sf_ini ! PUBLIC dom_vvl_sf ! !! * Module variables REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !: mut, muu, muv, muf !: REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time-step, = 2 rdttra ! ! except at nit000 (=rdttra) if neuler=0 !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Header$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS #if defined key_vvl SUBROUTINE dom_vvl_ini !!---------------------------------------------------------------------- !! *** ROUTINE dom_vvl_ini *** !! !! ** Purpose : compute coefficients muX at T-U-V-F points to spread !! ssh over the whole water column (scale factors) !! !!---------------------------------------------------------------------- INTEGER :: ji, jj, jk, zpk REAL(wp), DIMENSION(jpi,jpj) :: zmut, zmuu, zmuv, zmuf ! 2D workspace !!---------------------------------------------------------------------- IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'dom_vvl_ini : Variable volume activated' WRITE(numout,*) '~~~~~~~~~~~ compute coef. used to spread ssh over each layers' ENDIF IF( ln_zps ) CALL ctl_stop( 'dom_vvl_ini : option ln_zps is incompatible with variable volume option key_vvl') #if defined key_zco || defined key_dynspg_rl CALL ctl_stop( 'dom_vvl_ini : options key_zco/key_dynspg_rl are incompatible with variable volume option key_vvl') #endif fsvdept (:,:,:) = gdept (:,:,:) fsvdepw (:,:,:) = gdepw (:,:,:) fsvde3w (:,:,:) = gdep3w(:,:,:) fsve3t (:,:,:) = e3t (:,:,:) fsve3u (:,:,:) = e3u (:,:,:) fsve3v (:,:,:) = e3v (:,:,:) fsve3f (:,:,:) = e3f (:,:,:) fsve3w (:,:,:) = e3w (:,:,:) fsve3uw (:,:,:) = e3uw (:,:,:) fsve3vw (:,:,:) = e3vw (:,:,:) ! mu computation zmut(:,:) = 0.e0 zmuu(:,:) = 0.e0 zmuv(:,:) = 0.e0 zmuf(:,:) = 0.e0 mut (:,:,:) = 0.e0 muu (:,:,:) = 0.e0 muv (:,:,:) = 0.e0 muf (:,:,:) = 0.e0 DO jj = 1, jpj DO ji = 1, jpi zpk = mbathy(ji,jj) - 1 DO jk = 1, zpk zmut(ji,jj) = zmut(ji,jj) + fsve3t(ji,jj,jk) * SUM( fsve3t(ji,jj,jk:zpk) ) zmuu(ji,jj) = zmuu(ji,jj) + fsve3u(ji,jj,jk) * SUM( fsve3u(ji,jj,jk:zpk) ) zmuv(ji,jj) = zmuv(ji,jj) + fsve3v(ji,jj,jk) * SUM( fsve3v(ji,jj,jk:zpk) ) zmuf(ji,jj) = zmuf(ji,jj) + fsve3f(ji,jj,jk) * SUM( fsve3f(ji,jj,jk:zpk) ) END DO END DO END DO DO jj = 1, jpj DO ji = 1, jpi zpk = mbathy(ji,jj) - 1 DO jk = 1, zpk mut(ji,jj,jk) = SUM( fsve3t(ji,jj,jk:zpk) ) / zmut(ji,jj) muu(ji,jj,jk) = SUM( fsve3u(ji,jj,jk:zpk) ) / zmuu(ji,jj) muv(ji,jj,jk) = SUM( fsve3v(ji,jj,jk:zpk) ) / zmuv(ji,jj) muf(ji,jj,jk) = SUM( fsve3f(ji,jj,jk:zpk) ) / zmuf(ji,jj) END DO END DO END DO END SUBROUTINE dom_vvl_ini SUBROUTINE dom_vvl !!---------------------------------------------------------------------- !! *** ROUTINE dom_vvl *** !! !! ** Purpose : compute ssh at U-V-F points, T-W scale factors and local !! depths at each time step. !!---------------------------------------------------------------------- !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp), DIMENSION(jpi,jpj) :: zsshf ! 2D workspace !!---------------------------------------------------------------------- ! Sea Surface Height at u-v-fpoints DO jj = 1, jpjm1 DO ji = 1,jpim1 sshu(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) ! sshv(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) ! zsshf(ji,jj) = 0.25 * fmask(ji,jj,1) & & * ( sshn(ji ,jj) + sshn(ji ,jj+1) & & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) END DO END DO ! Boundaries conditions CALL lbc_lnk( sshu, 'U', 1. ) CALL lbc_lnk( sshv, 'V', 1. ) CALL lbc_lnk( zsshf, 'F', 1. ) ! Scale factors at T levels DO jk = 1, jpkm1 fse3t(:,:,jk) = fsve3t(:,:,jk) * ( 1 + sshn(:,:) * mut(:,:,jk) ) fse3u(:,:,jk) = fsve3u(:,:,jk) * ( 1 + sshu(:,:) * muu(:,:,jk) ) fse3v(:,:,jk) = fsve3v(:,:,jk) * ( 1 + sshv(:,:) * muv(:,:,jk) ) fse3f(:,:,jk) = fsve3f(:,:,jk) * ( 1 + zsshf(:,:) * muf(:,:,jk) ) END DO ! Scale factors at W levels fse3w (:,:,1) = fse3t(:,:,1) fse3uw(:,:,1) = fse3u(:,:,1) fse3vw(:,:,1) = fse3v(:,:,1) DO jk = 2, jpk fse3w (:,:,jk) = 0.5 * ( fse3t(:,:,jk-1) + fse3t(:,:,jk) ) fse3uw(:,:,jk) = 0.5 * ( fse3u(:,:,jk-1) + fse3u(:,:,jk) ) fse3vw(:,:,jk) = 0.5 * ( fse3v(:,:,jk-1) + fse3v(:,:,jk) ) END DO ! T and W points depth fsdept(:,:,1) = 0.5 * fse3w(:,:,1) fsdepw(:,:,1) = 0.e0 fsde3w(:,:,1) = fsdept(:,:,1) - sshn(:,:) DO jk = 2, jpk fsdept(:,:,jk) = fsdept(:,:,jk-1) + fse3w(:,:,jk) fsdepw(:,:,jk) = fsdepw(:,:,jk-1) + fse3t(:,:,jk-1) fsde3w(:,:,jk) = fsdept(:,:,jk ) - sshn(:,:) END DO ! Local depth or Inverse of the local depth of the water column at u- and v-points ! ------------------------------ ! Ocean depth at U- and V-points hu(:,:) = 0.e0 ; hv(:,:) = 0.e0 DO jk = 1, jpk hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk) hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) END DO ! Inverse of the local depth hur(:,:) = fse3u(:,:,1) ! Lower bound : thickness of the first model level hvr(:,:) = fse3v(:,:,1) DO jk = 2, jpk ! Sum of the vertical scale factors hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk) hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) END DO ! Compute and mask the inverse of the local depth hur(:,:) = 1. / hur(:,:) * umask(:,:,1) hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1) END SUBROUTINE dom_vvl SUBROUTINE dom_vvl_ssh( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dom_vvl_ssh *** !! !! ** Purpose : compute the ssha field used in tra_zdf, dynnxt, tranxt !! and dynspg routines !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! time step !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: z2dt, zraur ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zun, zvn, zhdiv ! 2D workspace !!---------------------------------------------------------------------- !! Sea surface height after (ssha) in variable volume case !! --------------------------====-----===============----- !! ssha is needed for the tracer time-stepping (trazdf_imp or !! tra_nxt), for the ssh time-stepping (dynspg_*) and for the !! dynamics time-stepping (dynspg_flt or dyn_nxt, and wzv). !! un, vn (or un_b and vn_b) and emp are needed, so it must be !! done after the call of oce_sbc. IF( neuler == 0 .AND. kt == nit000 ) THEN ! at nit000 r2dt(:) = rdttra(:) ! = rdtra (restarting with Euler time stepping) ELSEIF( kt <= nit000 + 1) THEN ! at nit000 or nit000+1 r2dt(:) = 2. * rdttra(:) ! = 2 rdttra (leapfrog) ENDIF z2dt = r2dt(1) zraur = 1. / rauw ! Vertically integrated quantities ! -------------------------------- IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN zun(:,:) = 0.e0 ; zvn(:,:) = 0.e0 ! DO jk = 1, jpkm1 ! ! Vertically integrated transports (now) zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) END DO ELSE ! lk_dynspg_ts is true zun (:,:) = un_b(:,:) ; zvn (:,:) = vn_b(:,:) ENDIF ! Horizontal divergence of barotropic transports !-------------------------------------------------- DO jj = 2, jpjm1 DO ji = 2, jpim1 ! vector opt. zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun(ji ,jj ) & & - e2u(ji-1,jj ) * zun(ji-1,jj ) & & + e1v(ji ,jj ) * zvn(ji ,jj ) & & - e1v(ji ,jj-1) * zvn(ji ,jj-1) ) & & / ( e1t(ji,jj) * e2t(ji,jj) ) END DO END DO #if defined key_obc && ( key_dynspg_exp || key_dynspg_ts ) ! open boundaries (div must be zero behind the open boundary) ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south #endif CALL lbc_lnk( zhdiv, 'T', 1. ) ! Sea surface elevation time stepping ! ----------------------------------- IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN ssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) ELSE ! lk_dynspg_ts is true ssha(:,:) = ( sshb_b(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) ENDIF END SUBROUTINE dom_vvl_ssh SUBROUTINE dom_vvl_sf( zssh, gridp, sfe3 ) !!---------------------------------------------------------------------- !! *** ROUTINE dom_vvl_sf *** !! !! ** Purpose : compute vertical scale factor at each time step !!---------------------------------------------------------------------- !! * Arguments CHARACTER(LEN=1) , INTENT( in ) :: gridp ! grid point type REAL(wp), DIMENSION(jpi,jpj) , INTENT( in ) :: zssh ! 2D workspace REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: sfe3 ! 3D workspace !! * Local declarations INTEGER :: jk ! dummy loop indices !!---------------------------------------------------------------------- SELECT CASE (gridp) CASE ('U') ! DO jk = 1, jpk sfe3(:,:,jk) = fsve3u(:,:,jk) * ( 1 + zssh(:,:) * muu(:,:,jk) ) ENDDO CASE ('V') ! DO jk = 1, jpk sfe3(:,:,jk) = fsve3v(:,:,jk) * ( 1 + zssh(:,:) * muv(:,:,jk) ) ENDDO CASE ('T') ! DO jk = 1, jpk sfe3(:,:,jk) = fsve3t(:,:,jk) * ( 1 + zssh(:,:) * mut(:,:,jk) ) ENDDO END SELECT END SUBROUTINE dom_vvl_sf SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini ) !!---------------------------------------------------------------------- !! *** ROUTINE dom_vvl_sf_ini *** !! !! ** Purpose : affect the appropriate vertical scale factor. It is done !! to avoid compiling error when using key_zco cpp key !!---------------------------------------------------------------------- !! * Arguments CHARACTER(LEN=1) , INTENT( in ) :: gridp ! grid point type REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini ! 3D workspace !!---------------------------------------------------------------------- SELECT CASE (gridp) CASE ('U') ! sfe3ini(:,:,:) = fse3u(:,:,:) CASE ('V') ! sfe3ini(:,:,:) = fse3v(:,:,:) CASE ('T') ! sfe3ini(:,:,:) = fse3t(:,:,:) END SELECT END SUBROUTINE dom_vvl_sf_ini #else !!---------------------------------------------------------------------- !! Default option : Empty routine !!---------------------------------------------------------------------- SUBROUTINE dom_vvl_ini END SUBROUTINE dom_vvl_ini SUBROUTINE dom_vvl END SUBROUTINE dom_vvl SUBROUTINE dom_vvl_ssh( kt ) !! * Arguments INTEGER, INTENT( in ) :: kt ! time step WRITE(*,*) 'dom_vvl_ssh: You should not have seen this print! error?', kt END SUBROUTINE dom_vvl_ssh SUBROUTINE dom_vvl_sf( zssh, gridp, sfe3 ) !! * Arguments CHARACTER(LEN=1) , INTENT( in ) :: gridp ! grid point type REAL(wp), DIMENSION(jpi,jpj) , INTENT( in ) :: zssh ! 2D workspace REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3 ! 3D workspace sfe3(:,:,:) = 0.e0 WRITE(*,*) 'sfe3: You should not have seen this print! error?', gridp WRITE(*,*) 'sfe3: You should not have seen this print! error?', zssh(1,1) END SUBROUTINE dom_vvl_sf SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini ) !! * Arguments CHARACTER(LEN=1) , INTENT( in ) :: gridp ! grid point type REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini ! 3D workspace sfe3ini(:,:,:) = 0.e0 WRITE(*,*) 'sfe3ini: You should not have seen this print! error?', gridp END SUBROUTINE dom_vvl_sf_ini #endif !!====================================================================== END MODULE domvvl