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 !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! 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 & e3u, & & e3v, & & e3f, & #endif & e1f, & & e2f, & & mig, & & mjg, & & nldi, & & nldj, & & nlei, & & nlej, & & umask, & & vmask USE dynadv , ONLY: & & ln_dynadv_vec ! vector form flag 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 ctl_stop ('vor_ene_tan not available yet') ! 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 ctl_stop ('vor_een_tan not available yet') ! CALL vor_een_tan( kt, ntot, ua_tl, va_tl ) ! total vorticity ! END SELECT END SUBROUTINE dyn_vor_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 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 ctl_stop ('vor_een_adj not available yet') ! 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, zuav, zvau ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 3D workspace 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_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 : oice 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 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) & & ) ! 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 ) !zun_tlin(:,:,:) = znu(:,:,:) !zvn_tlin(:,:,:) = znv(:,:,:) !zua_tlin(:,:,:) = zau(:,:,:) !zva_tlin(:,:,:) = zav(:,:,:) 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 ) !zrotn_tlin(:,:,:) = rotn_tl(:,:,:) 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(:,:,:) CALL dyn_vor_tan( nit000 ) 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(:,:,:) CALL dyn_vor_adj ( nitend ) 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' cl_name = 'dyn_vor_adj ' CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) 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