MODULE dynvor_tam #ifdef key_tam !!====================================================================== !! *** MODULE dynvor_tam *** !! Ocean dynamics: Update the momentum trend with the relative and !! planetary vorticity trends !! Tangent and Adjoint Module !!====================================================================== !! History of the drect module: !! 1.0 ! 89-12 (P. Andrich) vor_ens: Original code !! 5.0 ! 91-11 (G. Madec) vor_ene, vor_mix: Original code !! 6.0 ! 96-01 (G. Madec) s-coord, suppress work arrays !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !! 8.5 ! 04-02 (G. Madec) vor_een: Original code !! 9.0 ! 03-08 (G. Madec) vor_ctl: Original code !! 9.0 ! 05-11 (G. Madec) dyn_vor: Original code (new step architecture) !! 9.0 ! 06-11 (G. Madec) flux form advection: add metric term !! History of the TAM module: !! 9.0 ! 08-06 (A. Vidard) Skeleton !! 9.0 ! 09-01 (A. Vidard) TAM of the 06-11 version !! 9.0 ! 10-01 (F. Vigilant) Add een TAM option !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dyn_vor : Update the momentum trend with the vorticity trend !! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T) !! vor_ene : energy conserving scheme (ln_dynvor_ene=T) !! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T) !! vor_een : energy and enstrophy conserving (ln_dynvor_een=T) !! vor_ctl : set and control of the different vorticity option !!---------------------------------------------------------------------- USE par_kind, ONLY: & ! Precision variables & wp USE par_oce, ONLY: & ! Ocean space and time domain variables & jpi, & & jpj, & & jpk, & & jpim1, & & jpjm1, & & jpkm1, & & jpiglo USE oce , ONLY: & ! ocean dynamics and tracers & un, & & vn, & & rotn USE oce_tam , ONLY: & & un_tl, & & vn_tl, & & ua_tl, & & va_tl, & & rotn_tl, & & un_ad, & & vn_ad, & & ua_ad, & & va_ad, & & rotn_ad USE divcur , ONLY: & & div_cur USE divcur_tam , ONLY: & & div_cur_tan USE dom_oce , ONLY: & ! ocean space and time domain & ln_sco, & & ff, & & e1u, & & e2u, & & e1v, & & e2v, & #if defined key_zco & e3t_0, & #else & e3t, & & e3u, & & e3v, & & e3f, & #endif & e1f, & & e2f, & & mig, & & mjg, & & nldi, & & nldj, & & nlei, & & nlej, & & umask, & & vmask, & & tmask USE dynadv , ONLY: & & ln_dynadv_vec ! vector form flag USE lbclnk , ONLY: & ! Lateral boundary conditions & lbc_lnk USE in_out_manager, ONLY: & ! I/O manager & ctl_stop, & & lk_esopa, & & numnam, & & numout, & & nit000, & & nitend, & & lwp USE gridrandom , ONLY: & ! Random Gaussian noise on grids & grid_random USE dotprodfld, ONLY: & ! Computes dot product for 3D and 2D fields & dot_product USE tstool_tam , ONLY: & & prntst_adj, & ! ! random field standard deviation for: & stdu, & ! u-velocity & stdv ! v-velocity IMPLICIT NONE PRIVATE PUBLIC dyn_vor_tan ! routine called by step_tam.F90 PUBLIC dyn_vor_adj ! routine called by step_tam.F90 PUBLIC dyn_vor_adj_tst ! routine called by the tst.F90 !!* Namelist nam_dynvor: vorticity term LOGICAL, PUBLIC :: ln_dynvor_ene = .FALSE. !: energy conserving scheme LOGICAL, PUBLIC :: ln_dynvor_ens = .TRUE. !: enstrophy conserving scheme LOGICAL, PUBLIC :: ln_dynvor_mix = .FALSE. !: mixed scheme LOGICAL, PUBLIC :: ln_dynvor_een = .FALSE. !: energy and enstrophy conserving scheme INTEGER :: nvor = 0 ! type of vorticity trend used INTEGER :: ncor = 1 ! coriolis INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" CONTAINS SUBROUTINE dyn_vor_tan( kt ) !!---------------------------------------------------------------------- !! !! ** Purpose of the direct routine: !! compute the lateral ocean tracer physics. !! !! ** Action : - Update (ua,va) with the now vorticity term trend !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative !! and planetary vorticity trends) ('key_trddyn') !!---------------------------------------------------------------------- !! INTEGER, INTENT( in ) :: kt ! ocean time-step index !!---------------------------------------------------------------------- IF( kt == nit000 ) CALL vor_ctl_tam ! initialisation & control of options ! ! vorticity term SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend ! CASE ( -1 ) ! esopa: test all possibility with control print CALL vor_ene_tan( kt, ntot, ua_tl, va_tl ) CALL vor_ens_tan( kt, ntot, ua_tl, va_tl ) ! CALL vor_mix_tan( kt ) CALL vor_een_tan( kt, ntot, ua_tl, va_tl ) ! CASE ( 0 ) ! energy conserving scheme CALL vor_ene_tan( kt, ntot, ua_tl, va_tl ) ! total vorticity ! CASE ( 1 ) ! enstrophy conserving scheme CALL vor_ens_tan( kt, ntot, ua_tl, va_tl ) ! total vorticity ! CASE ( 2 ) ! mixed ene-ens scheme CALL ctl_stop ('vor_mix_tan not available yet') ! CALL vor_mix_tan( kt ) ! total vorticity (mix=ens-ene) ! CASE ( 3 ) ! energy and enstrophy conserving scheme CALL vor_een_tan( kt, ntot, ua_tl, va_tl ) ! total vorticity ! END SELECT END SUBROUTINE dyn_vor_tan SUBROUTINE vor_ene_tan( kt, kvor, pua_tl, pva_tl ) !!---------------------------------------------------------------------- !! *** ROUTINE vor_ene *** !! !! ** Purpose : Compute the now total vorticity trend and add it to !! the general trend of the momentum equation. !! !! ** Method : Trend evaluated using now fields (centered in time) !! and the Sadourny (1975) flux form formulation : conserves the !! horizontal kinetic energy. !! The trend of the vorticity term is given by: !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: !! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ] !! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ] !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: !! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ] !! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ] !! Add this trend to the general momentum trend (ua,va): !! (ua,va) = (ua,va) + ( voru , vorv ) !! !! ** Action : - Update (ua,va) with the now vorticity term trend !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative !! and planetary vorticity trends) ('key_trddyn') !! !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; ! ! =nrvm (relative vorticity or metric) REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_tl ! total u-trend REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_tl ! total v-trend !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zx1, zy1, zfact2 ! temporary scalars REAL(wp) :: zx2, zy2 ! " " REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 2D workspace REAL(wp) :: zx1tl, zy1tl ! temporary scalars REAL(wp) :: zx2tl, zy2tl ! " " REAL(wp), DIMENSION(jpi,jpj) :: zwxtl, zwytl, zwztl ! temporary 2D workspace !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF ! Local constant initialization zfact2 = 0.5 * 0.5 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) END DO END DO CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ff (ji,jj) & & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) END DO END DO END SELECT IF( ln_sco ) THEN zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) ELSE zwx(:,:) = e2u(:,:) * un(:,:,jk) zwy(:,:) = e1v(:,:) * vn(:,:,jk) ENDIF ! Tangent counterpart SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ; zwztl(:,:) = 0. ! planetary vorticity (Coriolis) CASE ( 2 ) ; zwztl(:,:) = rotn_tl(:,:,jk) ! relative vorticity CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwztl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) END DO END DO CASE ( 4 ) ; zwztl(:,:) = rotn_tl(:,:,jk) ! total (relative + planetary vorticity) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwztl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un_tl(ji ,jj+1,jk) + un_tl(ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) END DO END DO END SELECT IF( ln_sco ) THEN zwztl(:,:) = zwztl(:,:) / fse3f(:,:,jk) zwxtl(:,:) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk) zwytl(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk) ELSE zwxtl(:,:) = e2u(:,:) * un_tl(:,:,jk) zwytl(:,:) = e1v(:,:) * vn_tl(:,:,jk) ENDIF ! Compute and add the vorticity term trend ! ---------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) zy1tl = zwytl(ji,jj-1) + zwytl(ji+1,jj-1) zy2tl = zwytl(ji,jj ) + zwytl(ji+1,jj ) zx1tl = zwxtl(ji-1,jj) + zwxtl(ji-1,jj+1) zx2tl = zwxtl(ji ,jj) + zwxtl(ji ,jj+1) pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwztl(ji ,jj-1) * zy1 + zwz(ji ,jj-1) * zy1tl + zwztl(ji,jj) * zy2 + zwz(ji,jj) * zy2tl ) pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwztl(ji-1,jj ) * zx1 + zwz(ji-1,jj ) * zx1tl + zwztl(ji,jj) * zx2 + zwz(ji,jj) * zx2tl ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== END SUBROUTINE vor_ene_tan SUBROUTINE vor_ens_tan( kt, kvor, pua_tl, pva_tl ) !!---------------------------------------------------------------------- !! *** ROUTINE vor_ens_tan *** !! !! ** Purpose of the direct routine: !! Compute the now total vorticity trend and add it to !! the general trend of the momentum equation. !! !! ** Method of the direct routine: !! Trend evaluated using now fields (centered in time) !! and the Sadourny (1975) flux FORM formulation : conserves the !! potential enstrophy of a horizontally non-divergent flow. the !! trend of the vorticity term is given by: !! * s-coordinate (ln_sco=T), the e3. are inside the derivative: !! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] !! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ] !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: !! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ] !! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ] !! Add this trend to the general momentum trend (ua,va): !! (ua,va) = (ua,va) + ( voru , vorv ) !! !! ** Action : - Update (ua,va) arrays with the now vorticity term trend !! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative !! and planetary vorticity trends) ('key_trddyn') !! !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; ! ! =nrvm (relative vorticity or metric) REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_tl ! total u-trend REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_tl ! total v-trend !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zfact1, zuav, zvau ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 3D workspace REAL(wp) :: zuavtl, zvautl ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zwxtl, zwytl, zwztl ! temporary 3D workspace !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_vor_ens_tan : vorticity term: enstrophy conserving scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' ENDIF ! Local constant initialization zfact1 = 0.5 * 0.25 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) END DO END DO CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ff (ji,jj) & & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) END DO END DO END SELECT IF( ln_sco ) THEN DO jj = 1, jpj ! caution: don't use (:,:) for this loop DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) END DO END DO ELSE DO jj = 1, jpj ! caution: don't use (:,:) for this loop DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) END DO END DO ENDIF ! Compute and add the vorticity term trend ! ---------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== !CDIR PARALLEL DO PRIVATE( zwxtl, zwytl, zwztl ) ! =================== ! Tangent counterpart ! =================== ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ; zwztl(:,:) = 0.0_wp ! planetary vorticity (Coriolis) CASE ( 2 ,4) ; zwztl(:,:) = rotn_tl(:,:,jk) ! relative vorticity CASE ( 3 ,5 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwztl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) END DO END DO END SELECT IF( ln_sco ) THEN DO jj = 1, jpj ! caution: don't use (:,:) for this loop DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking zwztl(ji,jj) = zwztl(ji,jj) / fse3f(ji,jj,jk) zwxtl(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un_tl(ji,jj,jk) zwytl(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn_tl(ji,jj,jk) END DO END DO ELSE DO jj = 1, jpj ! caution: don't use (:,:) for this loop DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking zwxtl(ji,jj) = e2u(ji,jj) * un_tl(ji,jj,jk) zwytl(ji,jj) = e1v(ji,jj) * vn_tl(ji,jj,jk) END DO END DO ENDIF ! Compute and add the vorticity term trend ! ---------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zuavtl = zfact1 / e1u(ji,jj) * ( zwytl(ji ,jj-1) + zwytl(ji+1,jj-1) & & + zwytl(ji ,jj ) + zwytl(ji+1,jj ) ) zvautl =-zfact1 / e2v(ji,jj) * ( zwxtl(ji-1,jj ) + zwxtl(ji-1,jj+1) & & + zwxtl(ji ,jj ) + zwxtl(ji ,jj+1) ) pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk) & & + zuavtl * ( zwz( ji,jj-1) + zwz( ji,jj) ) & & + zuav * ( zwztl(ji,jj-1) + zwztl(ji,jj) ) pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk) & & + zvautl * ( zwz( ji-1,jj) + zwz( ji,jj) ) & & + zvau * ( zwztl(ji-1,jj) + zwztl(ji,jj) ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== END SUBROUTINE vor_ens_tan SUBROUTINE vor_een_tan( kt, kvor, pua_tl, pva_tl ) !!---------------------------------------------------------------------- !! *** ROUTINE vor_een_tan *** !! !! ** Purpose : Compute the now total vorticity trend and add it to !! the general trend of the momentum equation. !! !! ** Method : Trend evaluated using now fields (centered in time) !! and the Arakawa and Lamb (19XX) flux form formulation : conserves !! both the horizontal kinetic energy and the potential enstrophy !! when horizontal divergence is zero. !! The trend of the vorticity term is given by: !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: !! Add this trend to the general momentum trend (ua,va): !! (ua,va) = (ua,va) + ( voru , vorv ) !! !! ** Action : - Update (ua,va) with the now vorticity term trend !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative !! and planetary vorticity trends) ('key_trddyn') !! !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; ! ! =nrvm (relative vorticity or metric) REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_tl ! total u-trend REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_tl ! total v-trend !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zfac12, zua, zva ! temporary scalars REAL(wp) :: zuatl, zvatl ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 2D workspace REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse ! temporary 3D workspace REAL(wp), DIMENSION(jpi,jpj) :: zwxtl, zwytl, zwztl ! temporary 2D workspace REAL(wp), DIMENSION(jpi,jpj) :: ztnwtl, ztnetl, ztswtl, ztsetl ! temporary 3D workspace REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: ze3f !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn:vor_een_tam : vorticity term: energy and enstrophy conserving scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' DO jk = 1, jpk DO jj = 1, jpjm1 DO ji = 1, jpim1 ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) * 0.25_wp IF( ze3f(ji,jj,jk) /= 0.0_wp ) ze3f(ji,jj,jk) = 1.0_wp / ze3f(ji,jj,jk) END DO END DO END DO CALL lbc_lnk( ze3f, 'F', 1._wp ) ENDIF ! Local constant initialization zfac12 = 1.0_wp / 12.0_wp !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) zwz(:,:) = ff(:,:) * ze3f(:,:,jk) ! planetary vorticity (Coriolis) zwztl(:,:) = 0.0_wp CASE ( 2 ) zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) ! relative vorticity zwztl(:,:) = rotn_tl(:,:,jk) * ze3f(:,:,jk) CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) END DO END DO DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwztl(ji,jj) = ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) END DO END DO CASE ( 4 ) zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) ! total (relative + planetary vorticity) zwztl(:,:) = ( rotn_tl(:,:,jk) ) * ze3f(:,:,jk) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ff (ji,jj) & & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) * ze3f(ji,jj,jk) END DO END DO DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwztl(ji,jj) = ( ( ( vn_tl(ji+1,jj ,jk) + vn_tl (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un_tl(ji ,jj+1,jk) + un_tl (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) * ze3f(ji,jj,jk) END DO END DO END SELECT zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) zwxtl(:,:) = e2u(:,:) * fse3u(:,:,jk) * un_tl(:,:,jk) zwytl(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn_tl(:,:,jk) ! Compute and add the vorticity term trend ! ---------------------------------------- jj=2 ztne(1,:) = 0.0_wp ; ztnw(1,:) = 0.0_wp ; ztse(1,:) = 0.0_wp ; ztsw(1,:) = 0.0_wp ztnetl(1,:) = 0.0_wp ; ztnwtl(1,:) = 0.0_wp ; ztsetl(1,:) = 0.0_wp ; ztswtl(1,:) = 0.0_wp DO ji = 2, jpi ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ztnetl(ji,jj) = zwztl(ji-1,jj ) + zwztl(ji ,jj ) + zwztl(ji ,jj-1) ztnwtl(ji,jj) = zwztl(ji-1,jj-1) + zwztl(ji-1,jj ) + zwztl(ji ,jj ) ztsetl(ji,jj) = zwztl(ji ,jj ) + zwztl(ji ,jj-1) + zwztl(ji-1,jj-1) ztswtl(ji,jj) = zwztl(ji ,jj-1) + zwztl(ji-1,jj-1) + zwztl(ji-1,jj ) END DO DO jj = 3, jpj DO ji = fs_2, jpi ! vector opt. ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ztnetl(ji,jj) = zwztl(ji-1,jj ) + zwztl(ji ,jj ) + zwztl(ji ,jj-1) ztnwtl(ji,jj) = zwztl(ji-1,jj-1) + zwztl(ji-1,jj ) + zwztl(ji ,jj ) ztsetl(ji,jj) = zwztl(ji ,jj ) + zwztl(ji ,jj-1) + zwztl(ji-1,jj-1) ztswtl(ji,jj) = zwztl(ji ,jj-1) + zwztl(ji-1,jj-1) + zwztl(ji-1,jj ) END DO END DO DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zuatl = + zfac12 / e1u(ji,jj) * ( ztnetl(ji,jj ) * zwy(ji ,jj ) + ztne(ji,jj ) * zwytl(ji ,jj ) & & + ztnwtl(ji+1,jj) * zwy(ji+1,jj ) + ztnw(ji+1,jj) * zwytl(ji+1,jj ) & & + ztsetl(ji,jj ) * zwy(ji ,jj-1) + ztse(ji,jj ) * zwytl(ji ,jj-1) & & + ztswtl(ji+1,jj) * zwy(ji+1,jj-1) + ztsw(ji+1,jj) * zwytl(ji+1,jj-1)) zvatl = - zfac12 / e2v(ji,jj) * ( ztswtl(ji,jj+1) * zwx(ji-1,jj+1) + ztsw(ji,jj+1) * zwxtl(ji-1,jj+1) & & + ztsetl(ji,jj+1) * zwx(ji ,jj+1) + ztse(ji,jj+1) * zwxtl(ji ,jj+1) & & + ztnwtl(ji,jj ) * zwx(ji-1,jj ) + ztnw(ji,jj ) * zwxtl(ji-1,jj ) & & + ztnetl(ji,jj ) * zwx(ji ,jj ) + ztne(ji,jj ) * zwxtl(ji ,jj ) ) pua_tl(ji,jj,jk) = pua_tl(ji,jj,jk) + zuatl pva_tl(ji,jj,jk) = pva_tl(ji,jj,jk) + zvatl END DO END DO ! ! =============== END DO ! End of slab ! ! =============== END SUBROUTINE vor_een_tan SUBROUTINE dyn_vor_adj( kt ) !!---------------------------------------------------------------------- !! !! ** Purpose of the direct routine: !! compute the lateral ocean tracer physics. !! !! ** Action : - Update (ua,va) with the now vorticity term trend !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative !! and planetary vorticity trends) ('key_trddyn') !!---------------------------------------------------------------------- !! INTEGER, INTENT( in ) :: kt ! ocean time-step index !!---------------------------------------------------------------------- IF( kt == nitend ) CALL vor_ctl_tam ! initialisation & control of options ! ! vorticity term SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend ! CASE ( -1 ) ! esopa: test all possibility with control print CALL vor_een_adj( kt, ntot, ua_ad, va_ad ) ! CALL vor_mix_adj( kt ) CALL vor_ens_adj( kt, ntot, ua_ad, va_ad ) ! CALL vor_ene_adj( kt, ntot, ua_ad, va_ad ) ! CASE ( 0 ) ! energy conserving scheme CALL ctl_stop ('vor_ene_adj not available yet') ! CALL vor_ene_adj( kt, ntot, ua_ad, va_ad ) ! total vorticity ! CASE ( 1 ) ! enstrophy conserving scheme CALL vor_ens_adj( kt, ntot, ua_ad, va_ad ) ! total vorticity ! CASE ( 2 ) ! mixed ene-ens scheme CALL ctl_stop ('vor_mix_adj not available yet') ! CALL vor_mix_adj( kt ) ! total vorticity (mix=ens-ene) ! CASE ( 3 ) ! energy and enstrophy conserving scheme CALL vor_een_adj( kt, ntot, ua_ad, va_ad ) ! total vorticity ! END SELECT END SUBROUTINE dyn_vor_adj SUBROUTINE vor_ens_adj( kt, kvor, pua_ad, pva_ad ) !!---------------------------------------------------------------------- !! *** ROUTINE vor_ens_adj *** !! !! ** Purpose of the direct routine: !! Compute the now total vorticity trend and add it to !! the general trend of the momentum equation. !! !! ** Method of the direct routine: !! Trend evaluated using now fields (centered in time) !! and the Sadourny (1975) flux FORM formulation : conserves the !! potential enstrophy of a horizontally non-divergent flow. the !! trend of the vorticity term is given by: !! * s-coordinate (ln_sco=T), the e3. are inside the derivative: !! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] !! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ] !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: !! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ] !! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ] !! Add this trend to the general momentum trend (ua,va): !! (ua,va) = (ua,va) + ( voru , vorv ) !! !! ** Action : - Update (ua,va) arrays with the now vorticity term trend !! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative !! and planetary vorticity trends) ('key_trddyn') !! !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; ! ! =nrvm (relative vorticity or metric) REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_ad ! total u-trend REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_ad ! total v-trend !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zfact1 ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 3D workspace REAL(wp) :: zuav, zvau ! temporary scalars REAL(wp) :: zuavad, zvauad ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zwxad, zwyad, zwzad ! temporary 3D workspace !!---------------------------------------------------------------------- IF( kt == nitend ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_vor_ens_adj : vorticity term: enstrophy conserving scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' ENDIF ! Local constant initialization zfact1 = 0.5 * 0.25 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) )& & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) END DO END DO CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ff (ji,jj) & & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) END DO END DO END SELECT IF( ln_sco ) THEN DO jj = 1, jpj ! caution: don't use (:,:) for this loop DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) END DO END DO ELSE DO jj = 1, jpj ! caution: don't use (:,:) for this loop DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) END DO END DO ENDIF ! Compute and add the vorticity term trend ! ---------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== !CDIR PARALLEL DO PRIVATE( zwxad, zwyad, zwzad ) ! =================== ! Adjoint counterpart ! =================== zuavad = 0.0_wp zvauad = 0.0_wp zwxad(:,:) = 0.0_wp zwyad(:,:) = 0.0_wp zwzad(:,:) = 0.0_wp ! ! =============== DO jk = jpkm1, 1, -1 ! Horizontal slab ! ! =============== ! Compute and add the vorticity term trend ! ---------------------------------------- DO jj = jpjm1, 2, -1 DO ji = fs_jpim1, fs_2, -1 ! vector opt. zuavad = zuavad + pua_ad(ji,jj,jk) * ( zwz(ji,jj-1) + zwz(ji,jj) ) zwzad(ji,jj-1) = zwzad(ji,jj-1) + pua_ad(ji,jj,jk) * zuav zwzad(ji,jj ) = zwzad(ji,jj ) + pua_ad(ji,jj,jk) * zuav zvauad = zvauad + pva_ad(ji,jj,jk) * ( zwz(ji-1,jj) + zwz(ji,jj) ) zwzad(ji-1,jj) = zwzad(ji-1,jj) + pva_ad(ji,jj,jk) * zvau zwzad(ji ,jj) = zwzad(ji ,jj) + pva_ad(ji,jj,jk) * zvau zwyad(ji ,jj-1) = zwyad(ji ,jj-1) + zuavad * zfact1 / e1u(ji,jj) zwyad(ji+1,jj-1) = zwyad(ji+1,jj-1) + zuavad * zfact1 / e1u(ji,jj) zwyad(ji ,jj ) = zwyad(ji ,jj ) + zuavad * zfact1 / e1u(ji,jj) zwyad(ji+1,jj ) = zwyad(ji+1,jj ) + zuavad * zfact1 / e1u(ji,jj) zuavad = 0.0_wp zwxad(ji-1,jj ) = zwxad(ji-1,jj ) - zvauad * zfact1 / e2v(ji,jj) zwxad(ji-1,jj+1) = zwxad(ji-1,jj+1) - zvauad * zfact1 / e2v(ji,jj) zwxad(ji ,jj ) = zwxad(ji ,jj ) - zvauad * zfact1 / e2v(ji,jj) zwxad(ji ,jj+1) = zwxad(ji ,jj+1) - zvauad * zfact1 / e2v(ji,jj) zvauad = 0.0_wp END DO END DO IF( ln_sco ) THEN DO jj = jpj, 1, -1 ! caution: don't use (:,:) for this loop DO ji = jpi, 1, -1 ! it causes optimization problems on NEC in auto-tasking zwzad(ji,jj) = zwzad(ji,jj) / fse3f(ji,jj,jk) un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + zwxad(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + zwyad(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,jk) zwxad(ji,jj) = 0.0_wp zwyad(ji,jj) = 0.0_wp END DO END DO ELSE DO jj = jpj, 1, -1 ! caution: don't use (:,:) for this loop DO ji = jpi, 1, -1 ! it causes optimization problems on NEC in auto-tasking un_ad(ji,jj,jk) = un_ad(ji,jj,jk) + e2u(ji,jj) * zwxad(ji,jj) vn_ad(ji,jj,jk) = vn_ad(ji,jj,jk) + e1v(ji,jj) * zwyad(ji,jj) zwxad(ji,jj) = 0.0_wp zwyad(ji,jj) = 0.0_wp END DO END DO ENDIF ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ! planetary vorticity (Coriolis) zwzad(:,:) = 0.0_wp CASE ( 2 ,4) ! relative vorticity rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:) zwzad(:,:) = 0.0_wp CASE ( 3 ,5 ) ! metric term DO jj = jpjm1, 1, -1 DO ji = fs_jpim1, 1, -1 ! vector opt. vn_ad(ji+1,jj,jk) = vn_ad(ji+1,jj,jk) & & + zwzad(ji,jj) * ( e2v(ji+1,jj) - e2v(ji,jj) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) vn_ad(ji ,jj,jk) = vn_ad(ji ,jj,jk) & & + zwzad(ji,jj) * ( e2v(ji+1,jj) - e2v(ji,jj) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) un_ad(ji,jj+1,jk) = un_ad(ji,jj+1,jk) & & - zwzad(ji,jj) * ( e1u(ji,jj+1) - e1u(ji,jj) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) un_ad(ji,jj ,jk) = un_ad(ji,jj ,jk) & & - zwzad(ji,jj) * ( e1u(ji,jj+1) - e1u(ji,jj) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) zwzad(ji,jj) = 0.0_wp END DO END DO END SELECT ! ! =============== END DO ! End of slab ! ! =============== END SUBROUTINE vor_ens_adj SUBROUTINE vor_een_adj( kt, kvor, pua_ad, pva_ad ) !!---------------------------------------------------------------------- !! *** ROUTINE vor_een_adj *** !! !! ** Purpose : Compute the now total vorticity trend and add it to !! the general trend of the momentum equation. !! !! ** Method : Trend evaluated using now fields (centered in time) !! and the Arakawa and Lamb (19XX) flux form formulation : conserves !! both the horizontal kinetic energy and the potential enstrophy !! when horizontal divergence is zero. !! The trend of the vorticity term is given by: !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: !! Add this trend to the general momentum trend (ua,va): !! (ua,va) = (ua,va) + ( voru , vorv ) !! !! ** Action : - Update (ua,va) with the now vorticity term trend !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative !! and planetary vorticity trends) ('key_trddyn') !! !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; ! ! =nrvm (relative vorticity or metric) REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua_ad ! total u-trend REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva_ad ! total v-trend !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zfac12 ! temporary scalars REAL(wp) :: zuaad, zvaad ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 2D workspace REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse ! temporary 3D workspace REAL(wp), DIMENSION(jpi,jpj) :: zwxad, zwyad, zwzad ! temporary 2D workspace REAL(wp), DIMENSION(jpi,jpj) :: ztnwad, ztnead, ztswad, ztsead ! temporary 3D workspace REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: ze3f !!---------------------------------------------------------------------- ! local adjoint initailization zuaad = 0.0_wp ; zvaad = 0.0_wp zwxad (:,:) = 0.0_wp ; zwyad (:,:) = 0.0_wp ; zwzad (:,:) = 0.0_wp ztnwad(:,:) = 0.0_wp ; ztnead(:,:) = 0.0_wp ; ztswad(:,:) = 0.0_wp ; ztsead(:,:) = 0.0_wp IF( kt == nitend ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn:vor_een_adj : vorticity term: energy and enstrophy conserving scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' DO jk = 1, jpk DO jj = 1, jpjm1 DO ji = 1, jpim1 ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) * 0.25_wp IF( ze3f(ji,jj,jk) /= 0.0_wp ) ze3f(ji,jj,jk) = 1.0_wp / ze3f(ji,jj,jk) END DO END DO END DO CALL lbc_lnk( ze3f, 'F', 1._wp ) ENDIF ! Local constant initialization zfac12 = 1.0_wp / 12.0_wp !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Potential vorticity and horizontal fluxes (Direct local variables init) ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ; zwz(:,:) = ff(:,:) * ze3f(:,:,jk) ! planetary vorticity (Coriolis) CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) ! relative vorticity CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) )& & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) END DO END DO CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) ! total (relative + planetary vorticity) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ff (ji,jj) & & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) * ze3f(ji,jj,jk) END DO END DO END SELECT zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) ! Compute and add the vorticity term trend ! ---------------------------------------- jj=2 ztne(1,:) = 0.0_wp ; ztnw(1,:) = 0.0_wp ; ztse(1,:) = 0.0_wp ; ztsw(1,:) = 0.0_wp DO ji = 2, jpi ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) END DO DO jj = 3, jpj DO ji = fs_2, jpi ! vector opt. ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) END DO END DO ! =================== ! Adjoint counterpart ! =================== DO jj = jpjm1, 2, -1 DO ji = fs_jpim1, fs_2, -1 ! vector opt. zuaad = zuaad + pua_ad(ji,jj,jk) zvaad = zvaad + pva_ad(ji,jj,jk) zvaad = - zvaad * zfac12 / e2v(ji,jj) ztswad(ji ,jj+1) = ztswad(ji ,jj+1) + zvaad * zwx (ji-1,jj+1) zwxad (ji-1,jj+1) = zwxad (ji-1,jj+1) + zvaad * ztsw(ji ,jj+1) ztsead(ji ,jj+1) = ztsead(ji ,jj+1) + zvaad * zwx (ji ,jj+1) zwxad (ji ,jj+1) = zwxad (ji ,jj+1) + zvaad * ztse(ji ,jj+1) ztnwad(ji ,jj ) = ztnwad(ji ,jj ) + zvaad * zwx (ji-1,jj ) zwxad (ji-1,jj ) = zwxad (ji-1,jj ) + zvaad * ztnw(ji ,jj ) ztnead(ji ,jj ) = ztnead(ji ,jj ) + zvaad * zwx (ji ,jj ) zwxad (ji ,jj ) = zwxad (ji ,jj ) + zvaad * ztne(ji ,jj ) zvaad = 0.0_wp zuaad = zuaad * zfac12 / e1u(ji,jj) ztnead(ji ,jj ) = ztnead(ji ,jj ) + zuaad * zwy (ji ,jj ) zwyad (ji ,jj ) = zwyad (ji ,jj ) + zuaad * ztne(ji ,jj ) ztnwad(ji+1,jj ) = ztnwad(ji+1,jj ) + zuaad * zwy (ji+1,jj ) zwyad (ji+1,jj ) = zwyad (ji+1,jj ) + zuaad * ztnw(ji+1,jj ) ztsead(ji ,jj ) = ztsead(ji ,jj ) + zuaad * zwy (ji ,jj-1) zwyad (ji ,jj-1) = zwyad (ji ,jj-1) + zuaad * ztse(ji ,jj ) ztswad(ji+1,jj ) = ztswad(ji+1,jj ) + zuaad * zwy (ji+1,jj-1) zwyad (ji+1,jj-1) = zwyad (ji+1,jj-1) + zuaad * ztsw(ji+1,jj ) zuaad = 0.0_wp END DO END DO DO jj = jpj, 3, -1 DO ji = jpi, fs_2, -1 ! vector opt. zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztswad(ji,jj) zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztswad(ji,jj) zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztswad(ji,jj) ztswad(ji ,jj ) = 0.0_wp zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztsead(ji,jj) zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztsead(ji,jj) zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztsead(ji,jj) ztsead(ji,jj) = 0.0_wp zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztnwad(ji,jj) zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztnwad(ji,jj) zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztnwad(ji,jj) ztnwad(ji ,jj ) = 0.0_wp zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztnead(ji,jj) zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztnead(ji,jj) zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztnead(ji,jj) ztnead(ji,jj) = 0.0_wp END DO END DO jj=2 DO ji = jpi, 2, -1 zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztswad(ji,jj) zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztswad(ji,jj) zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztswad(ji,jj) ztswad(ji,jj) = 0.0_wp zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztsead(ji,jj) zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztsead(ji,jj) zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztsead(ji,jj) ztsead(ji ,jj ) = 0.0_wp zwzad (ji-1,jj-1) = zwzad(ji-1,jj-1) + ztnwad(ji,jj) zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztnwad(ji,jj) zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztnwad(ji,jj) ztnwad(ji ,jj ) = 0.0_wp zwzad (ji-1,jj ) = zwzad(ji-1,jj ) + ztnead(ji,jj) zwzad (ji ,jj ) = zwzad(ji ,jj ) + ztnead(ji,jj) zwzad (ji ,jj-1) = zwzad(ji ,jj-1) + ztnead(ji,jj) ztnead(ji ,jj ) = 0.0_wp END DO ztnead(1,:) = 0.0_wp ; ztnwad(1,:) = 0.0_wp ztsead(1,:) = 0.0_wp ; ztswad(1,:) = 0.0_wp vn_ad(:,:,jk) = vn_ad(:,:,jk) + zwyad(:,:) * e1v(:,:) * fse3v(:,:,jk) un_ad(:,:,jk) = un_ad(:,:,jk) + zwxad(:,:) * e2u(:,:) * fse3u(:,:,jk) zwyad(:,:) = 0.0_wp zwxad(:,:) = 0.0_wp ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) zwzad(:,:) = 0.0_wp CASE ( 2 ) rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:) * ze3f(:,:,jk) zwzad(:,:) = 0.0_wp CASE ( 3 ) ! metric term DO jj = jpjm1, 1, -1 DO ji = fs_jpim1, 1, -1 ! vector opt. zwzad(ji ,jj ) = zwzad(ji,jj) * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) vn_ad(ji+1,jj ,jk) = vn_ad(ji+1,jj ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) vn_ad(ji ,jj ,jk) = vn_ad(ji ,jj ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) un_ad(ji ,jj+1,jk) = - un_ad(ji ,jj+1,jk) + zwzad(ji,jj) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) un_ad(ji ,jj ,jk) = - un_ad(ji ,jj ,jk) + zwzad(ji,jj) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) zwzad(ji ,jj ) = 0.0_wp END DO END DO CASE ( 4 ) rotn_ad(:,:,jk) = rotn_ad(:,:,jk) + zwzad(:,:) * ze3f(:,:,jk) zwzad(:,:) = 0.0_wp CASE ( 5 ) ! total (coriolis + metric) DO jj = jpjm1, 1, -1 DO ji = fs_jpim1, 1, -1 ! vector opt. zwzad(ji ,jj ) = zwzad(ji,jj) * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) vn_ad(ji+1,jj ,jk) = vn_ad(ji+1,jj ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) vn_tl(ji ,jj ,jk) = vn_tl(ji ,jj ,jk) + zwzad(ji,jj) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) un_ad(ji ,jj+1,jk) = un_ad(ji ,jj+1,jk) - zwzad(ji,jj) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) un_ad(ji ,jj ,jk) = un_ad(ji ,jj ,jk) - zwzad(ji,jj) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) zwzad(ji ,jj ) = 0.0_wp END DO END DO END SELECT ! ! =============== END DO ! End of slab ! ! =============== END SUBROUTINE vor_een_adj SUBROUTINE vor_ctl_tam !!--------------------------------------------------------------------- !! *** ROUTINE vor_ctl_tam *** !! !! ** Purpose : Control the consistency between cpp options for !! tracer advection schemes !!---------------------------------------------------------------------- INTEGER :: ioptio ! temporary integer NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een !!---------------------------------------------------------------------- REWIND ( numnam ) ! Read Namelist nam_dynvor : Vorticity scheme options READ ( numnam, nam_dynvor ) IF(lwp) THEN ! Namelist print WRITE(numout,*) WRITE(numout,*) 'dyn:vor_ctl_tam : vorticity term : read namelist and control the consistency' WRITE(numout,*) '~~~~~~~~~~~~~~~' WRITE(numout,*) ' Namelist nam_dynvor : choice of the vorticity term scheme' WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een ENDIF ioptio = 0 ! Control of vorticity scheme options IF( ln_dynvor_ene ) ioptio = ioptio + 1 IF( ln_dynvor_ens ) ioptio = ioptio + 1 IF( ln_dynvor_mix ) ioptio = ioptio + 1 IF( ln_dynvor_een ) ioptio = ioptio + 1 IF( lk_esopa ) ioptio = 1 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) ! ! Set nvor (type of scheme for vorticity) IF( ln_dynvor_ene ) nvor = 0 IF( ln_dynvor_ens ) nvor = 1 IF( ln_dynvor_mix ) nvor = 2 IF( ln_dynvor_een ) nvor = 3 IF( lk_esopa ) nvor = -1 ! ! Set ncor, nrvm, ntot (type of vorticity) IF(lwp) WRITE(numout,*) ncor = 1 IF( ln_dynadv_vec ) THEN IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity' nrvm = 2 ntot = 4 ELSE IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term' nrvm = 3 ntot = 5 ENDIF IF(lwp) THEN ! Print the choice WRITE(numout,*) IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme' IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme' IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme' IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme' IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options' ENDIF ! END SUBROUTINE vor_ctl_tam SUBROUTINE dyn_vor_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE dyn_adv_adj_tst *** !! !! ** Purpose : Test the adjoint routine. !! !! ** Method : Verify the scalar product !! !! ( L dx )^T W dy = dx^T L^T W dy !! !! where L = tangent routine !! L^T = adjoint routine !! W = diagonal matrix of scale factors !! dx = input perturbation (random field) !! dy = L dx !! !! ** Action : Separate tests are applied for the following dx and dy: !! !! 1) dx = ( SSH ) and dy = ( SSH ) !! !! History : !! ! 08-08 (A. Vidard) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit INTEGER :: & & ji, & ! dummy loop indices & jj, & & jk, & & jt INTEGER, DIMENSION(jpi,jpj) :: & & iseed_2d ! 2D seed for the random number generator !! * Local declarations REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zun_tlin, & ! Tangent input: now u-velocity & zvn_tlin, & ! Tangent input: now v-velocity & zrotn_tlin, & ! Tangent input: now rot & zun_adout, & ! Adjoint output: now u-velocity & zvn_adout, & ! Adjoint output: now v-velocity & zrotn_adout, & ! Adjoint output: now rot & zua_adout, & ! Tangent output: after u-velocity & zva_adout, & ! Tangent output: after v-velocity & zua_tlin, & ! Tangent output: after u-velocity & zva_tlin, & ! Tangent output: after v-velocity & zua_tlout, & ! Tangent output: after u-velocity & zva_tlout, & ! Tangent output: after v-velocity & zua_adin, & ! Tangent output: after u-velocity & zva_adin, & ! Tangent output: after v-velocity & zau, & ! 3D random field for rotn & zav, & ! 3D random field for rotn & znu, & ! 3D random field for u & znv ! 3D random field for v REAL(KIND=wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp1_1, & ! scalar product components & zsp1_2, & & zsp2, & ! scalar product involving the adjoint routine & zsp2_1, & ! scalar product components & zsp2_2, & & zsp2_3, & & zsp2_4, & & zsp2_5 CHARACTER(LEN=14) :: cl_name ! Allocate memory ALLOCATE( & & zun_tlin(jpi,jpj,jpk), & & zvn_tlin(jpi,jpj,jpk), & & zrotn_tlin(jpi,jpj,jpk), & & zun_adout(jpi,jpj,jpk), & & zvn_adout(jpi,jpj,jpk), & & zrotn_adout(jpi,jpj,jpk), & & zua_adout(jpi,jpj,jpk), & & zva_adout(jpi,jpj,jpk), & & zua_tlin(jpi,jpj,jpk), & & zva_tlin(jpi,jpj,jpk), & & zua_tlout(jpi,jpj,jpk), & & zva_tlout(jpi,jpj,jpk), & & zua_adin(jpi,jpj,jpk), & & zva_adin(jpi,jpj,jpk), & & zau(jpi,jpj,jpk), & & zav(jpi,jpj,jpk), & & znu(jpi,jpj,jpk), & & znv(jpi,jpj,jpk) & & ) ! init ntot parameter CALL vor_ctl_tam ! initialisation & control of options DO jt = 1, 2 IF (jt == 1) nvor=1 ! enstrophy conserving scheme IF (jt == 2) nvor=3 ! energy and enstrophy conserving scheme ! Initialize rotn CALL div_cur ( nit000 ) !================================================================== ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and ! dy = ( hdivb_tl, hdivn_tl ) !================================================================== !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zun_tlin(:,:,:) = 0.0_wp zvn_tlin(:,:,:) = 0.0_wp zrotn_tlin(:,:,:) = 0.0_wp zun_adout(:,:,:) = 0.0_wp zvn_adout(:,:,:) = 0.0_wp zrotn_adout(:,:,:) = 0.0_wp zua_tlout(:,:,:) = 0.0_wp zva_tlout(:,:,:) = 0.0_wp zua_adin(:,:,:) = 0.0_wp zva_adin(:,:,:) = 0.0_wp zua_adout(:,:,:) = 0.0_wp zva_adout(:,:,:) = 0.0_wp zua_tlin(:,:,:) = 0.0_wp zva_tlin(:,:,:) = 0.0_wp znu(:,:,:) = 0.0_wp znv(:,:,:) = 0.0_wp zau(:,:,:) = 0.0_wp zav(:,:,:) = 0.0_wp un_tl(:,:,:) = 0.0_wp vn_tl(:,:,:) = 0.0_wp ua_tl(:,:,:) = 0.0_wp va_tl(:,:,:) = 0.0_wp un_ad(:,:,:) = 0.0_wp vn_ad(:,:,:) = 0.0_wp ua_ad(:,:,:) = 0.0_wp va_ad(:,:,:) = 0.0_wp rotn_tl(:,:,:) = 0.0_wp rotn_ad(:,:,:) = 0.0_wp !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 596035 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znu, 'U', 0.0_wp, stdu ) DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 523432 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znv, 'V', 0.0_wp, stdv ) DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 432545 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, zau, 'U', 0.0_wp, stdu ) DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 287503 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, zav, 'V', 0.0_wp, stdv ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zun_tlin(ji,jj,jk) = znu(ji,jj,jk) zvn_tlin(ji,jj,jk) = znv(ji,jj,jk) zua_tlin(ji,jj,jk) = zau(ji,jj,jk) zva_tlin(ji,jj,jk) = zav(ji,jj,jk) END DO END DO END DO un_tl(:,:,:) = zun_tlin(:,:,:) vn_tl(:,:,:) = zvn_tlin(:,:,:) ua_tl(:,:,:) = zua_tlin(:,:,:) va_tl(:,:,:) = zva_tlin(:,:,:) ! initialize rotn_tl with noise CALL div_cur_tan ( nit000 ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zrotn_tlin(ji,jj,jk) = rotn_tl(ji,jj,jk) END DO END DO END DO rotn_tl(:,:,:) = zrotn_tlin(:,:,:) IF (nvor == 1 ) CALL vor_ens_tan( nit000, ntot, ua_tl, va_tl ) IF (nvor == 3 ) CALL vor_een_tan( nit000, ntot, ua_tl, va_tl ) zua_tlout(:,:,:) = ua_tl(:,:,:) zva_tlout(:,:,:) = va_tl(:,:,:) !-------------------------------------------------------------------- ! Initialize the adjoint variables: dy^* = W dy !-------------------------------------------------------------------- DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) & & * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) & & * umask(ji,jj,jk) zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) & & * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) & & * vmask(ji,jj,jk) END DO END DO END DO !-------------------------------------------------------------------- ! Compute the scalar product: ( L dx )^T W dy !-------------------------------------------------------------------- zsp1_1 = DOT_PRODUCT( zua_tlout, zua_adin ) zsp1_2 = DOT_PRODUCT( zva_tlout, zva_adin ) zsp1 = zsp1_1 + zsp1_2 !-------------------------------------------------------------------- ! Call the adjoint routine: dx^* = L^T dy^* !-------------------------------------------------------------------- ua_ad(:,:,:) = zua_adin(:,:,:) va_ad(:,:,:) = zva_adin(:,:,:) IF (nvor == 1 ) CALL vor_ens_adj( nitend, ntot, ua_ad, va_ad ) IF (nvor == 3 ) CALL vor_een_adj( nitend, ntot, ua_ad, va_ad ) zun_adout(:,:,:) = un_ad(:,:,:) zvn_adout(:,:,:) = vn_ad(:,:,:) zrotn_adout(:,:,:) = rotn_ad(:,:,:) zua_adout(:,:,:) = ua_ad(:,:,:) zva_adout(:,:,:) = va_ad(:,:,:) zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout ) zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout ) zsp2_3 = DOT_PRODUCT( zrotn_tlin, zrotn_adout ) zsp2_4 = DOT_PRODUCT( zua_tlin, zua_adout ) zsp2_5 = DOT_PRODUCT( zva_tlin, zva_adout ) zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 ! Compare the scalar products ! 14 char:'12345678901234' IF (nvor == 1 ) cl_name = 'dynvor_adj ens' IF (nvor == 3 ) cl_name = 'dynvor_adj een' CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) END DO DEALLOCATE( & & zun_tlin, & & zvn_tlin, & & zrotn_tlin, & & zun_adout, & & zvn_adout, & & zrotn_adout, & & zua_adout, & & zva_adout, & & zua_tlin, & & zva_tlin, & & zua_tlout, & & zva_tlout, & & zua_adin, & & zva_adin, & & zau, & & zav, & & znu, & & znv & & ) END SUBROUTINE dyn_vor_adj_tst !!============================================================================= #endif END MODULE dynvor_tam