MODULE dynspg_flt_tam !!---------------------------------------------------------------------- !! This software is governed by the CeCILL licence (Version 2) !!---------------------------------------------------------------------- #if defined key_tam !!====================================================================== !! *** MODULE dynspg_flt_tam : TANGENT/ADJOINT OF MODULE dynspg_flt *** !! !! Ocean dynamics: surface pressure gradient trend !! !!====================================================================== !! History of the direct module: !! History OPA ! 1998-05 (G. Roullet) free surface !! ! 1998-10 (G. Madec, M. Imbard) release 8.2 !! NEMO O.1 ! 2002-08 (G. Madec) F90: Free form and module !! - ! 2002-11 (C. Talandier, A-M Treguier) Open boundaries !! 1.0 ! 2004-08 (C. Talandier) New trends organization !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization !! 2.0 ! 2006-07 (S. Masson) distributed restart using iom !! - ! 2006-08 (J.Chanut, A.Sellar) Calls to BDY routines. !! 3.2 ! 2009-03 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module !! History of the TAM module: !! 9.0 ! 2008-08 (A. Vidard) skeleton !! - ! 2008-08 (A. Weaver) original version !! NEMO 3.2 ! 2010-04 (F. Vigilant) converison to 3.2 !!---------------------------------------------------------------------- # if defined key_dynspg_flt || defined key_esopa !!---------------------------------------------------------------------- !! 'key_dynspg_flt' filtered free surface !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dyn_spg_flt_tan : update the momentum trend with the surface pressure !! gradient in the filtered free surface case !! (tangent routine) !! dyn_spg_flt_adj : update the momentum trend with the surface pressure !! gradient in the filtered free surface case !! (adjoint routine) !! dyn_spg_flt_adj_tst : Test of the adjoint routine !!---------------------------------------------------------------------- USE par_kind , ONLY: & ! Precision variables & wp USE prtctl ! Print control USE par_oce , ONLY: & ! Ocean space and time domain variables & jpi, & & jpj, & & jpk, & & jpim1, & & jpjm1, & & jpkm1, & & jpr2di, & & jpr2dj, & & jpij, & & jpiglo, & & lk_vopt_loop USE in_out_manager, ONLY: & ! I/O manager & lwp, & & numout, & & nit000, & & nitend, & & ctl_stop, & & ln_ctl, & & ctmp1 USE phycst , ONLY: & ! Physical constants & grav USE lib_mpp , ONLY: & ! distributed memory computing & lk_mpp, & & mpp_sum USE mppsum , ONLY: & ! Summation of arrays across processors & mpp_sum_indep USE dom_oce , ONLY: & ! Ocean space and time domain & e1u, & & e2u, & & e1v, & & e2v, & & e1t, & & e2t, & # if defined key_zco & e3t_0, & # else & e3u, & & e3v, & # endif & tmask, & & umask, & & vmask, & & bmask, & & rdt, & & neuler, & & n_cla, & & mig, & & mjg, & & nldi, & & nldj, & & nlei, & & nlej, & & atfp, & & atfp1, & & lk_vvl USE solver , ONLY: & ! Solvers & sol_mat USE sol_oce , ONLY: & ! Solver variables in memory & nn_solv, & & nn_nmod, & & ncut, & & niter, & & eps, & & epsr, & & rnorme, & & res, & & c_solver_pt, & & gcdprc USE oce_tam , ONLY: & ! Tangent-linear and adjoint model variables & ub_tl, & & vb_tl, & & ua_tl, & & va_tl, & & ub_ad, & & vb_ad, & & ua_ad, & & va_ad, & & spgu_tl, & & spgv_tl, & & spgu_ad, & & spgv_ad, & & sshn_tl, & & sshn_ad USE sbc_oce_tam , ONLY: & ! Surface BCs for tangent and adjoint model & emp_tl, & & emp_ad # if defined key_obc # if defined key_pomme_r025 USE obc_oce USE obcdyn_tam # else Error, OBC not ready. # endif # endif USE sol_oce , ONLY: & ! Solver arrays & gcdmat USE sol_oce_tam , ONLY: & ! Solver arrays for tangent and adjoint model & gcx_tl, & & gcxb_tl, & & gcb_tl, & & gcx_ad, & & gcxb_ad, & & gcb_ad USE solsor_tam , ONLY: & ! SOR solver for tangent and adjoint model & sol_sor_tan, & & sol_sor_adj USE solpcg_tam , ONLY: & ! PCG solver for tangent and adjoint model & sol_pcg_tan, & & sol_pcg_adj USE lbclnk , ONLY: & ! Lateral boundary conditions & lbc_lnk, & & lbc_lnk_e USE lbclnk_tam , ONLY: & ! TAM lateral boundary conditions & lbc_lnk_adj, & & lbc_lnk_e_adj USE gridrandom , ONLY: & ! Random Gaussian noise on grids & grid_random USE dotprodfld , ONLY: & ! Computes dot product for 3D and 2D fields & dot_product USE paresp , ONLY: & ! Normalized energy weights & wesp_emp, & & wesp_ssh USE dynadv, ONLY: & & ln_dynadv_vec ! vector form flag USE cla_dynspg_tam, ONLY: & ! Cross land advection on surface pressure & dyn_spg_cla_tan, & & dyn_spg_cla_adj USE tstool_tam , ONLY: & & stdu, & & stdv, & & stdw, & & stdemp, & & stdssh, & & prntst_adj, & & prntst_tlm IMPLICIT NONE PRIVATE !! * Accessibility PUBLIC dyn_spg_flt_tan, & ! routine called by step_tan.F90 & dyn_spg_flt_adj, & ! routine called by step_adj.F90 & dyn_spg_flt_adj_tst ! routine called by the tst.F90 # if defined key_tst_tlm PUBLIC dyn_spg_flt_tlm_tst # endif !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_spg_flt_tan( kt, kindic ) !!--------------------------------------------------------------------- !! *** routine dyn_spg_flt_tan *** !! !! ** Purpose of the direct routine: !! Compute the now trend due to the surface pressure !! gradient in case of filtered free surface formulation and add !! it to the general trend of momentum equation. !! !! ** Method of the direct routine: !! Filtered free surface formulation. The surface !! pressure gradient is given by: !! spgu = 1/rau0 d/dx(ps) = 1/e1u di( sshn + btda ) !! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + btda ) !! where sshn is the free surface elevation and btda is the after !! time derivative of the free surface elevation !! -1- evaluate the surface presure trend (including the addi- !! tional force) in three steps: !! a- compute the right hand side of the elliptic equation: !! gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ] !! where (spgu,spgv) are given by: !! spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ] !! - grav 2 rdt hu /e1u di[sshn + emp] !! spgv = vertical sum[ e3v (vb+ 2 rdt va) ] !! - grav 2 rdt hv /e2v dj[sshn + emp] !! and define the first guess from previous computation : !! zbtd = btda !! btda = 2 zbtd - btdb !! btdb = zbtd !! b- compute the relative accuracy to be reached by the !! iterative solver !! c- apply the solver by a call to sol... routine !! -2- compute and add the free surface pressure gradient inclu- !! ding the additional force used to stabilize the equation. !! !! ** Action : - Update (ua,va) with the surf. pressure gradient trend !! !! References : Roullet and Madec 1999, JGR. !!--------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( IN ) :: & & kt ! ocean time-step index INTEGER, INTENT( OUT ) :: & & kindic ! solver convergence flag (<0 if not converge) !! * Local declarations INTEGER :: & & ji, & ! dummy loop indices & jj, & & jk REAL(wp) :: & & z2dt, & ! temporary scalars & z2dtg, & & zgcb, & & zbtd, & & ztdgu, & & ztdgv !! Variable volume REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace & zsshub, & & zsshua, & & zsshvb, & & zsshva, & & zssha REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D workspace & zfse3ub, & & zfse3ua, & & zfse3vb, & & zfse3va !!---------------------------------------------------------------------- ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_spg_flt_tan : surface pressure gradient trend' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~ (free surface constant volume case)' ! set to zero free surface specific arrays spgu_tl(:,:) = 0.0_wp ! surface pressure gradient (i-direction) spgv_tl(:,:) = 0.0_wp ! surface pressure gradient (j-direction) ENDIF ! Local constant initialization IF ( neuler == 0 .AND. kt == nit000 ) THEN z2dt = rdt ! time step: Euler if restart from rest CALL sol_mat(kt) ! initialize matrix ELSEIF ( neuler == 0 .AND. kt == nit000 + 1 ) THEN z2dt = 2.0_wp * rdt ! time step: leap-frog CALL sol_mat(kt) ! initialize matrix ELSEIF ( neuler /= 0 .AND. kt == nit000 ) THEN z2dt = 2.0_wp * rdt ! time step: leap-frog CALL sol_mat(kt) ! initialize matrix ELSE z2dt = 2.0_wp * rdt ! time step: leap-frog ENDIF z2dtg = grav * z2dt ! Evaluate the masked next velocity (effect of the additional force not included) IF( lk_vvl ) THEN ! variable volume (surface pressure gradient already included in dyn_hpg) ! IF( ln_dynadv_vec ) THEN ! vector form : applied on velocity DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua_tl(ji,jj,jk) = ( ub_tl(ji,jj,jk) + z2dt * ua_tl(ji,jj,jk) ) * umask(ji,jj,jk) va_tl(ji,jj,jk) = ( vb_tl(ji,jj,jk) + z2dt * va_tl(ji,jj,jk) ) * vmask(ji,jj,jk) END DO END DO END DO ! ELSE ! flux form : applied on thickness weighted velocity DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua_tl(ji,jj,jk) = ( ub_tl(ji,jj,jk) * fse3u_b(ji,jj,jk) & & + z2dt * ua_tl(ji,jj,jk) * fse3u_n(ji,jj,jk) ) & & / fse3u_a(ji,jj,jk) * umask(ji,jj,jk) va_tl(ji,jj,jk) = ( vb_tl(ji,jj,jk) * fse3v_b(ji,jj,jk) & & + z2dt * va_tl(ji,jj,jk) * fse3v_n(ji,jj,jk) ) & & / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk) END DO END DO END DO ! ENDIF ! ELSE ! fixed volume (add the surface pressure gradient + unweighted time stepping) ! DO jj = 2, jpjm1 ! Surface pressure gradient (now) DO ji = fs_2, fs_jpim1 ! vector opt. spgu_tl(ji,jj) = - grav * ( sshn_tl(ji+1,jj) - sshn_tl(ji,jj) ) / e1u(ji,jj) spgv_tl(ji,jj) = - grav * ( sshn_tl(ji,jj+1) - sshn_tl(ji,jj) ) / e2v(ji,jj) END DO END DO DO jk = 1, jpkm1 ! unweighted time stepping DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua_tl(ji,jj,jk) = ( ub_tl(ji,jj,jk) + z2dt * ( ua_tl(ji,jj,jk) + spgu_tl(ji,jj) ) ) * umask(ji,jj,jk) va_tl(ji,jj,jk) = ( vb_tl(ji,jj,jk) + z2dt * ( va_tl(ji,jj,jk) + spgv_tl(ji,jj) ) ) * vmask(ji,jj,jk) END DO END DO END DO ! ENDIF # if defined key_obc IF( lk_obc ) CALL obc_dyn_tan( kt ) ! Update velocities on each open boundary with the radiation algorithm ! IF( lk_obc ) CALL obc_vol_tan( kt ) ! Correction of the barotropic componant velocity to control the volume of the system # endif # if defined key_bdy CALL ctl_stop( 'dyn_spg_flt_tan: key_bdy is not available' ) # endif # if defined key_agrif CALL ctl_stop( 'dyn_spg_flt_tan: key_agrif is not available' ) # endif # if defined key_orca_r2 IF( n_cla == 1 ) CALL dyn_spg_cla_tan( kt ) ! Cross Land Advection (update (ua,va)) # endif ! compute the next vertically averaged velocity (effect of the additional force not included) ! --------------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. spgu_tl(ji,jj) = 0.0_wp spgv_tl(ji,jj) = 0.0_wp END DO END DO ! vertical sum !CDIR NOLOOPCHG IF( lk_vopt_loop ) THEN ! vector opt., forced unroll DO jk = 1, jpkm1 DO ji = 1, jpij spgu_tl(ji,1) = spgu_tl(ji,1) + fse3u(ji,1,jk) * ua_tl(ji,1,jk) spgv_tl(ji,1) = spgv_tl(ji,1) + fse3v(ji,1,jk) * va_tl(ji,1,jk) END DO END DO ELSE ! No vector opt. DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 spgu_tl(ji,jj) = spgu_tl(ji,jj) + fse3u(ji,jj,jk) * ua_tl(ji,jj,jk) spgv_tl(ji,jj) = spgv_tl(ji,jj) + fse3v(ji,jj,jk) * va_tl(ji,jj,jk) END DO END DO END DO ENDIF ! transport: multiplied by the horizontal scale factor DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. spgu_tl(ji,jj) = spgu_tl(ji,jj) * e2u(ji,jj) spgv_tl(ji,jj) = spgv_tl(ji,jj) * e1v(ji,jj) END DO END DO CALL lbc_lnk( spgu_tl, 'U', -1.0_wp ) ! lateral boundary conditions CALL lbc_lnk( spgv_tl, 'V', -1.0_wp ) IF( lk_vvl ) CALL ctl_stop( 'dyn_spg_flt_tan: lk_vvl is not available' ) ! Right hand side of the elliptic equation and first guess ! ----------------------------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! Divergence of the after vertically averaged velocity zgcb = spgu_tl(ji,jj) - spgu_tl(ji-1,jj) & & + spgv_tl(ji,jj) - spgv_tl(ji,jj-1) gcb_tl(ji,jj) = gcdprc(ji,jj) * zgcb ! First guess of the after barotropic transport divergence zbtd = gcx_tl(ji,jj) gcx_tl (ji,jj) = 2.0_wp * zbtd - gcxb_tl(ji,jj) gcxb_tl(ji,jj) = zbtd END DO END DO ! apply the lateral boundary conditions IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) CALL lbc_lnk_e( gcb_tl, c_solver_pt, 1.0_wp ) # if defined key_agrif CALL ctl_stop( 'dyn_spg_flt_tan: key_agrif is not available' ) # endif ! Relative precision (computation on one processor) ! ------------------ rnorme = SUM( gcb_tl(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb_tl(1:jpi,1:jpj) * bmask(:,:) ) IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain epsr = eps * eps * rnorme ncut = 0 ! if rnorme is 0, the solution is 0, the solver is not called IF( rnorme == 0.0_wp ) THEN gcx_tl(:,:) = 0.0_wp res = 0.0_wp niter = 0 ncut = 999 ENDIF ! Evaluate the next transport divergence ! -------------------------------------- ! Iterarive solver for the elliptic equation (except IF sol.=0) ! (output in gcx with boundary conditions applied) kindic = 0 IF( ncut == 0 ) THEN IF ( nn_solv == 1 ) THEN ; CALL sol_pcg_tan( kt, kindic ) ! diagonal preconditioned conjuguate gradient ELSEIF( nn_solv == 2 ) THEN ; CALL sol_sor_tan( kt, kindic ) ! successive-over-relaxation ENDIF ENDIF ! Transport divergence gradient multiplied by z2dt ! --------------------------------------------==== DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! trend of Transport divergence gradient ztdgu = z2dtg * ( gcx_tl(ji+1,jj ) - gcx_tl(ji,jj) ) / e1u(ji,jj) ztdgv = z2dtg * ( gcx_tl(ji ,jj+1) - gcx_tl(ji,jj) ) / e2v(ji,jj) ! multiplied by z2dt # if defined key_obc ! caution : grad D = 0 along open boundaries IF( Agrif_Root() ) THEN spgu_tl(ji,jj) = z2dt * ztdgu * obcumask(ji,jj) spgv_tl(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj) ELSE spgu_tl(ji,jj) = z2dt * ztdgu spgv_tl(ji,jj) = z2dt * ztdgv ENDIF # elif defined key_bdy CALL ctl_stop( 'dyn_spg_flt_tan: key_bdy is not available' ) # else spgu_tl(ji,jj) = z2dt * ztdgu spgv_tl(ji,jj) = z2dt * ztdgv # endif END DO END DO # if defined key_agrif CALL ctl_stop( 'dyn_spg_flt_tan: key_agrif is not available' ) # endif ! Add the trends multiplied by z2dt to the after velocity ! ------------------------------------------------------- ! ( c a u t i o n : (ua_tl,va_tl) here are the after velocity not the ! trend, the leap-frog time stepping will not ! be done in dynnxt_tan.F90 routine) DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua_tl(ji,jj,jk) = ( ua_tl(ji,jj,jk) + spgu_tl(ji,jj) ) * umask(ji,jj,jk) va_tl(ji,jj,jk) = ( va_tl(ji,jj,jk) + spgv_tl(ji,jj) ) * vmask(ji,jj,jk) END DO END DO END DO ! END SUBROUTINE dyn_spg_flt_tan SUBROUTINE dyn_spg_flt_adj( kt, kindic ) !!---------------------------------------------------------------------- !! *** routine dyn_spg_flt_adj *** !! !! ** Purpose of the direct routine: !! Compute the now trend due to the surface pressure !! gradient in case of filtered free surface formulation and add !! it to the general trend of momentum equation. !! !! ** Method of the direct routine: !! Filtered free surface formulation. The surface !! pressure gradient is given by: !! spgu = 1/rau0 d/dx(ps) = 1/e1u di( sshn + btda ) !! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + btda ) !! where sshn is the free surface elevation and btda is the after !! time derivative of the free surface elevation !! -1- evaluate the surface presure trend (including the addi- !! tional force) in three steps: !! a- compute the right hand side of the elliptic equation: !! gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ] !! where (spgu,spgv) are given by: !! spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ] !! - grav 2 rdt hu /e1u di[sshn + emp] !! spgv = vertical sum[ e3v (vb+ 2 rdt va) ] !! - grav 2 rdt hv /e2v dj[sshn + emp] !! and define the first guess from previous computation : !! zbtd = btda !! btda = 2 zbtd - btdb !! btdb = zbtd !! b- compute the relative accuracy to be reached by the !! iterative solver !! c- apply the solver by a call to sol... routine !! -2- compute and add the free surface pressure gradient inclu- !! ding the additional force used to stabilize the equation. !! !! ** Action : - Update (ua,va) with the surf. pressure gradient trend !! !! References : Roullet and Madec 1999, JGR. !!--------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( IN ) :: & & kt ! ocean time-step index INTEGER, INTENT( OUT ) :: & & kindic ! solver convergence flag (<0 if not converge) !! * Local declarations INTEGER :: & & ji, & ! dummy loop indices & jj, & & jk REAL(wp) :: & & z2dt, & ! temporary scalars & z2dtg, & & zgcb, & & zbtd, & & ztdgu, & & ztdgv !!---------------------------------------------------------------------- ! IF( kt == nitend ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_spg_flt_adj : surface pressure gradient trend' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~ (free surface constant volume case)' ENDIF ! Local constant initialization IF ( neuler == 0 .AND. kt == nit000 ) THEN z2dt = rdt ! time step: Euler if restart from rest CALL sol_mat(kt) ! initialize matrix ELSEIF ( neuler == 0 .AND. kt == nit000 + 1 ) THEN z2dt = 2.0_wp * rdt ! time step: leap-frog CALL sol_mat(kt) ! reinitialize matrix ELSEIF ( kt == nitend ) THEN z2dt = 2.0_wp * rdt ! time step: leap-frog CALL sol_mat(kt) ! reinitialize matrix ELSEIF ( neuler /= 0 .AND. kt == nit000 ) THEN z2dt = 2.0_wp * rdt ! time step: leap-frog CALL sol_mat(kt) ! initialize matrix ELSE z2dt = 2.0_wp * rdt ! time step: leap-frog ENDIF z2dtg = grav * z2dt ! Add the trends multiplied by z2dt to the after velocity ! ----------------------------------------------------------- ! ( c a u t i o n : (ua_ad,va_ad) here are the after velocity not the ! trend, the leap-frog time stepping will not ! be done in dynnxt_adj.F90 routine) DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) * umask(ji,jj,jk) va_ad(ji,jj,jk) = va_ad(ji,jj,jk) * vmask(ji,jj,jk) spgu_ad(ji,jj) = spgu_ad(ji,jj) + ua_ad(ji,jj,jk) spgv_ad(ji,jj) = spgv_ad(ji,jj) + va_ad(ji,jj,jk) END DO END DO END DO # if defined key_agrif CALL ctl_stop( 'dyn_spg_flt_adj: key_agrif is not available' ) # endif ! Transport divergence gradient multiplied by z2dt ! --------------------------------------------==== DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! multiplied by z2dt # if defined key_obc ! caution : grad D = 0 along open boundaries IF( Agrif_Root() ) THEN ztdgu = z2dt * spgu_ad(ji,jj) * obcumask(ji,jj) ztdgv = z2dt * spgv_ad(ji,jj) * obcvmask(ji,jj) spgu_ad(ji,jj) = 0.0_wp spgv_ad(ji,jj) = 0.0_wp ELSE ztdgu = z2dt * spgu_ad(ji,jj) ztdgv = z2dt * spgv_ad(ji,jj) spgu_ad(ji,jj) = 0.0_wp spgv_ad(ji,jj) = 0.0_wp ENDIF # elif defined key_bdy CALL ctl_stop( 'dyn_spg_flt_adj: key_bdy is not available' ) # else ztdgu = z2dt * spgu_ad(ji,jj) ztdgv = z2dt * spgv_ad(ji,jj) spgu_ad(ji,jj) = 0.0_wp spgv_ad(ji,jj) = 0.0_wp # endif ! trend of Transport divergence gradient ztdgu = ztdgu * z2dtg / e1u(ji,jj) ztdgv = ztdgv * z2dtg / e2v(ji,jj) gcx_ad(ji ,jj ) = gcx_ad(ji ,jj ) - ztdgu - ztdgv gcx_ad(ji ,jj+1) = gcx_ad(ji ,jj+1) + ztdgv gcx_ad(ji+1,jj ) = gcx_ad(ji+1,jj ) + ztdgu END DO END DO ! Evaluate the next transport divergence ! -------------------------------------- ! Iterative solver for the elliptic equation (except IF sol.=0) ! (output in gcx_ad with boundary conditions applied) kindic = 0 ncut = 0 ! Force solver IF( ncut == 0 ) THEN IF ( nn_solv == 1 ) THEN ; CALL sol_pcg_adj( kt, kindic ) ! diagonal preconditioned conjuguate gradient ELSEIF( nn_solv == 2 ) THEN ; CALL sol_sor_adj( kt, kindic ) ! successive-over-relaxation ENDIF ENDIF # if defined key_agrif CALL ctl_stop( 'dyn_spg_flt_adj: key_agrif is not available' ) # endif ! Right hand side of the elliptic equation and first guess ! -------------------------------------------------------- ! apply the lateral boundary conditions IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) & & CALL lbc_lnk_e_adj( gcb_ad, c_solver_pt, 1.0_wp ) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! First guess of the after barotropic transport divergence zbtd = gcxb_ad(ji,jj) + 2.0_wp * gcx_ad(ji,jj) gcxb_ad(ji,jj) = - gcx_ad(ji,jj) gcx_ad (ji,jj) = zbtd ! Divergence of the after vertically averaged velocity zgcb = gcb_ad(ji,jj) * gcdprc(ji,jj) gcb_ad(ji,jj) = 0.0_wp spgu_ad(ji-1,jj ) = spgu_ad(ji-1,jj ) - zgcb spgu_ad(ji ,jj ) = spgu_ad(ji ,jj ) + zgcb spgv_ad(ji ,jj-1) = spgv_ad(ji ,jj-1) - zgcb spgv_ad(ji ,jj ) = spgv_ad(ji ,jj ) + zgcb END DO END DO IF( lk_vvl ) CALL ctl_stop( 'dyn_spg_flt_adj: lk_vvl is not available' ) ! Boundary conditions on (spgu_ad,spgv_ad) CALL lbc_lnk_adj( spgu_ad, 'U', -1.0_wp ) CALL lbc_lnk_adj( spgv_ad, 'V', -1.0_wp ) ! transport: multiplied by the horizontal scale factor DO jj = 2,jpjm1 DO ji = fs_2,fs_jpim1 ! vector opt. spgu_ad(ji,jj) = spgu_ad(ji,jj) * e2u(ji,jj) spgv_ad(ji,jj) = spgv_ad(ji,jj) * e1v(ji,jj) END DO END DO ! compute the next vertically averaged velocity (effect of the additional force not included) ! --------------------------------------------- ! vertical sum !CDIR NOLOOPCHG IF( lk_vopt_loop ) THEN ! vector opt., forced unroll DO jk = 1, jpkm1 DO ji = 1, jpij ua_ad(ji,1,jk) = ua_ad(ji,1,jk) + fse3u(ji,1,jk) * spgu_ad(ji,1) va_ad(ji,1,jk) = va_ad(ji,1,jk) + fse3v(ji,1,jk) * spgv_ad(ji,1) END DO END DO ELSE ! No vector opt. DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = 2, jpim1 ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) + fse3u(ji,jj,jk) * spgu_ad(ji,jj) va_ad(ji,jj,jk) = va_ad(ji,jj,jk) + fse3v(ji,jj,jk) * spgv_ad(ji,jj) END DO END DO END DO ENDIF DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. spgu_ad(ji,jj) = 0.0_wp spgv_ad(ji,jj) = 0.0_wp END DO END DO # if defined key_orca_r2 IF( n_cla == 1 ) CALL dyn_spg_cla_adj( kt ) ! Cross Land Advection (update (ua,va)) # endif # if defined key_agrif CALL ctl_stop( 'dyn_spg_flt_adj: key_agrif is not available' ) # endif # if defined key_bdy CALL ctl_stop( 'dyn_spg_flt_adj: key_bdy is not available' ) # endif # if defined key_obc ! IF( lk_obc ) CALL obc_vol_adj( kt ) ! Correction of the barotropic componant velocity to control the volume of the system IF( lk_obc ) CALL obc_dyn_adj( kt ) ! Update velocities on each open boundary with the radiation algorithm # endif ! Evaluate the masked next velocity (effect of the additional force not included) IF( lk_vvl ) THEN ! variable volume (surface pressure gradient already included in dyn_hpg) ! IF( ln_dynadv_vec ) THEN ! vector form : applied on velocity DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ub_ad(ji,jj,jk) = ub_ad(ji,jj,jk) + z2dt * ua_ad(ji,jj,jk) * umask(ji,jj,jk) ua_ad(ji,jj,jk) = z2dt * ua_ad(ji,jj,jk) * umask(ji,jj,jk) vb_ad(ji,jj,jk) = vb_ad(ji,jj,jk) + z2dt * va_ad(ji,jj,jk) * vmask(ji,jj,jk) va_ad(ji,jj,jk) = z2dt * va_ad(ji,jj,jk) * vmask(ji,jj,jk) END DO END DO END DO ! ELSE ! flux form : applied on thickness weighted velocity DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) / fse3u_a(ji,jj,jk) * umask(ji,jj,jk) ub_ad(ji,jj,jk) = ub_ad(ji,jj,jk) + ua_ad(ji,jj,jk) * fse3u_b(ji,jj,jk) ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) * z2dt * fse3u_n(ji,jj,jk) va_ad(ji,jj,jk) = va_ad(ji,jj,jk) / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk) vb_ad(ji,jj,jk) = vb_ad(ji,jj,jk) + va_ad(ji,jj,jk) * fse3v_b(ji,jj,jk) va_ad(ji,jj,jk) = va_ad(ji,jj,jk) * z2dt * fse3v_n(ji,jj,jk) END DO END DO END DO ! ENDIF ! ELSE ! fixed volume (add the surface pressure gradient + unweighted time stepping) ! DO jk = 1, jpkm1 ! unweighted time stepping DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua_ad( ji,jj,jk) = ua_ad(ji,jj,jk) * umask(ji,jj,jk) ub_ad( ji,jj,jk) = ub_ad(ji,jj,jk) + ua_ad(ji,jj,jk) spgu_ad(ji,jj ) = spgu_ad(ji,jj) + ua_ad(ji,jj,jk) * z2dt ua_ad( ji,jj,jk) = ua_ad(ji,jj,jk) * z2dt va_ad( ji,jj,jk) = va_ad(ji,jj,jk) * vmask(ji,jj,jk) vb_ad( ji,jj,jk) = vb_ad(ji,jj,jk) + va_ad(ji,jj,jk) spgv_ad(ji,jj ) = spgv_ad(ji,jj) + va_ad(ji,jj,jk) * z2dt va_ad( ji,jj,jk) = va_ad(ji,jj,jk) * z2dt END DO END DO END DO DO jj = 2, jpjm1 ! Surface pressure gradient (now) DO ji = fs_2, fs_jpim1 ! vector opt. spgu_ad(ji ,jj ) = spgu_ad(ji ,jj ) * grav / e1u(ji,jj) spgv_ad(ji ,jj ) = spgv_ad(ji ,jj ) * grav / e2v(ji,jj) sshn_ad(ji ,jj ) = sshn_ad(ji ,jj ) + spgv_ad(ji,jj) sshn_ad(ji ,jj+1) = sshn_ad(ji ,jj+1) - spgv_ad(ji,jj) sshn_ad(ji ,jj ) = sshn_ad(ji ,jj ) + spgu_ad(ji,jj) sshn_ad(ji+1,jj ) = sshn_ad(ji+1,jj ) - spgu_ad(ji,jj) spgu_ad(ji ,jj ) = 0.0_wp spgv_ad(ji ,jj ) = 0.0_wp END DO END DO ENDIF IF( kt == nit000 ) THEN ! set to zero free surface specific arrays spgu_ad(:,:) = 0.0_wp ! surface pressure gradient (i-direction) spgv_ad(:,:) = 0.0_wp ! surface pressure gradient (j-direction) ENDIF END SUBROUTINE dyn_spg_flt_adj SUBROUTINE dyn_spg_flt_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE dyn_spg_flt_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 : !! !! History : !! ! 09-01 (A. Weaver) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit !! * Local declarations REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zua_tlin, & ! Tangent input: ua_tl & zva_tlin, & ! Tangent input: va_tl & zub_tlin, & ! Tangent input: ub_tl & zvb_tlin, & ! Tangent input: vb_tl & zua_tlout, & ! Tangent output: ua_tl & zva_tlout, & ! Tangent output: va_tl & zub_tlout, & ! Tangent output: ua_tl & zvb_tlout, & ! Tangent output: va_tl & zub_adin, & ! Tangent output: ua_ad & zvb_adin, & ! Tangent output: va_ad & zua_adin, & ! Adjoint input: ua_ad & zva_adin, & ! Adjoint input: va_ad & zua_adout, & ! Adjoint output: ua_ad & zva_adout, & ! Adjoint output: va_ad & zub_adout, & ! Adjoint oputput: ub_ad & zvb_adout, & ! Adjoint output: vb_ad & znu ! 3D random field for u REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & zgcx_tlin, zgcxb_tlin, zgcb_tlin, zgcx_tlout, zgcxb_tlout, zgcb_tlout, & & zgcx_adin, zgcxb_adin, zgcb_adin, zgcx_adout, zgcxb_adout, zgcb_adout, & & zspgu_tlout, zspgv_tlout, zspgu_adin, zspgv_adin REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & & zsshn_tlin, & ! Tangent input: sshn_tl & zsshn_adout,& ! Adjoint output: sshn_ad & zemp_tlin, & ! Tangent input: emp_tl & zemp_adout, & ! Adjoint output: emp_ad & znssh ! 2D random field for SSH REAL(wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp2 ! scalar product involving the adjoint routine INTEGER, DIMENSION(jpi,jpj) :: & & iseed_2d ! 2D seed for the random number generator INTEGER :: & & indic, & & istp INTEGER :: & & ji, & & jj, & & jk, & & kmod, & & jstp CHARACTER (LEN=14) :: & & cl_name INTEGER :: & & jpert INTEGER, PARAMETER :: jpertmax = 7 ! Allocate memory ALLOCATE( & & zua_tlin(jpi,jpj,jpk), & & zva_tlin(jpi,jpj,jpk), & & zub_tlin(jpi,jpj,jpk), & & zvb_tlin(jpi,jpj,jpk), & & zua_tlout(jpi,jpj,jpk), & & zva_tlout(jpi,jpj,jpk), & & zua_adin(jpi,jpj,jpk), & & zva_adin(jpi,jpj,jpk), & & zua_adout(jpi,jpj,jpk), & & zva_adout(jpi,jpj,jpk), & & zub_adout(jpi,jpj,jpk), & & zvb_adout(jpi,jpj,jpk), & & znu(jpi,jpj,jpk) & & ) ALLOCATE( & & zsshn_tlin(jpi,jpj), & & zsshn_adout(jpi,jpj),& & zemp_tlin(jpi,jpj), & & zemp_adout(jpi,jpj), & & znssh(jpi,jpj) & & ) # if defined key_obc ALLOCATE( zemp_tlout(jpi,jpj), zemp_adin (jpi,jpj), & zwn_tlout (jpi,jpj), zwn_adin (jpi,jpj) ) # endif ALLOCATE( zgcx_tlin (jpi,jpj), zgcx_tlout (jpi,jpj), zgcx_adin (jpi,jpj), zgcx_adout (jpi,jpj), & zgcxb_tlin(jpi,jpj), zgcxb_tlout(jpi,jpj), zgcxb_adin(jpi,jpj), zgcxb_adout(jpi,jpj), & zgcb_tlin (jpi,jpj), zgcb_tlout (jpi,jpj), zgcb_adin (jpi,jpj), zgcb_adout (jpi,jpj) & & ) ALLOCATE ( zub_tlout(jpi,jpj,jpk), zvb_tlout(jpi,jpj,jpk), & zub_adin (jpi,jpj,jpk), zvb_adin (jpi,jpj,jpk) ) ALLOCATE( zspgu_tlout (jpi,jpj), zspgv_tlout (jpi,jpj), zspgu_adin (jpi,jpj), zspgv_adin (jpi,jpj)) !========================================================================= ! dx = ( sshb_tl, sshn_tl, ub_tl, ua_tl, vb_tl, va_tl, wn_tl, emp_tl ) ! and dy = ( sshb_tl, sshn_tl, ua_tl, va_tl ) !========================================================================= ! Test for time steps nit000 and nit000 + 1 (the matrix changes) DO jstp = nit000, nit000 + 1 DO jpert = 1, jpertmax istp = jstp !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zub_tlin (:,:,:) = 0.0_wp zvb_tlin (:,:,:) = 0.0_wp zua_tlin (:,:,:) = 0.0_wp zva_tlin (:,:,:) = 0.0_wp zua_tlout(:,:,:) = 0.0_wp zva_tlout(:,:,:) = 0.0_wp zua_adin (:,:,:) = 0.0_wp zva_adin (:,:,:) = 0.0_wp zub_adout(:,:,:) = 0.0_wp zvb_adout(:,:,:) = 0.0_wp zua_adout(:,:,:) = 0.0_wp zva_adout(:,:,:) = 0.0_wp zsshn_tlin (:,:) = 0.0_wp zemp_tlin (:,:) = 0.0_wp zsshn_adout(:,:) = 0.0_wp zemp_adout (:,:) = 0.0_wp zspgu_adin (:,:) = 0.0_wp zspgv_adin (:,:) = 0.0_wp zspgu_tlout(:,:) = 0.0_wp zspgv_tlout(:,:) = 0.0_wp zgcx_tlout (:,:) = 0.0_wp ; zgcx_adin (:,:) = 0.0_wp ; zgcx_adout (:,:) = 0.0_wp zgcxb_tlout(:,:) = 0.0_wp ; zgcxb_adin(:,:) = 0.0_wp ; zgcxb_adout(:,:) = 0.0_wp zgcb_tlout (:,:) = 0.0_wp ; zgcb_adin (:,:) = 0.0_wp ; zgcb_adout (:,:) = 0.0_wp ub_tl(:,:,:) = 0.0_wp vb_tl(:,:,:) = 0.0_wp ua_tl(:,:,:) = 0.0_wp va_tl(:,:,:) = 0.0_wp sshn_tl(:,:) = 0.0_wp emp_tl(:,:) = 0.0_wp gcb_tl(:,:) = 0.0_wp gcx_tl(:,:) = 0.0_wp gcxb_tl(:,:) = 0.0_wp spgu_tl(:,:) = 0.0_wp spgv_tl(:,:) = 0.0_wp ub_ad(:,:,:) = 0.0_wp vb_ad(:,:,:) = 0.0_wp ua_ad(:,:,:) = 0.0_wp va_ad(:,:,:) = 0.0_wp sshn_ad(:,:) = 0.0_wp emp_ad(:,:) = 0.0_wp gcb_ad(:,:) = 0.0_wp gcx_ad(:,:) = 0.0_wp gcxb_ad(:,:) = 0.0_wp spgu_ad(:,:) = 0.0_wp spgv_ad(:,:) = 0.0_wp !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- IF ( (jpert == 1) .OR. (jpert == jpertmax) ) THEN DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 456953 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znu, 'U', 0.0_wp, stdu ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zua_tlin(ji,jj,jk) = znu(ji,jj,jk) END DO END DO END DO ENDIF IF ( (jpert == 2) .OR. (jpert == jpertmax) ) THEN DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 3434334 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znu, 'V', 0.0_wp, stdv ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zva_tlin(ji,jj,jk) = znu(ji,jj,jk) END DO END DO END DO ENDIF IF ( (jpert == 3) .OR. (jpert == jpertmax) ) THEN DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 935678 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znu, 'U', 0.0_wp, stdu ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zub_tlin(ji,jj,jk) = znu(ji,jj,jk) END DO END DO END DO ENDIF IF ( (jpert == 4) .OR. (jpert == jpertmax) ) THEN DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 1462578 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znu, 'V', 0.0_wp, stdv ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zvb_tlin(ji,jj,jk) = znu(ji,jj,jk) END DO END DO END DO ENDIF IF ( (jpert == 5) .OR. (jpert == jpertmax) ) THEN DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 25674910 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdemp ) DO jj = nldj, nlej DO ji = nldi, nlei zemp_tlin(ji,jj) = znssh(ji,jj) END DO END DO ENDIF IF ( (jpert == 6) .OR. (jpert == jpertmax) ) THEN DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 12672456 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdssh ) DO jj = nldj, nlej DO ji = nldi, nlei zsshn_tlin(ji,jj) = znssh(ji,jj) END DO END DO END IF # if defined key_obc zgcx_tlin (:,:) = ( zua_tlin(:,:,1) + zub_tlin(:,:,1) ) / 10. zgcxb_tlin (:,:) = ( zua_tlin(:,:,2) + zub_tlin(:,:,2) ) / 10. zgcb_tlin (:,:) = ( zua_tlin(:,:,3) + zub_tlin(:,:,3) ) / 10. # endif zgcx_tlin (:,:) = ( zua_tlin(:,:,1) + zub_tlin(:,:,1) ) / 10. zgcxb_tlin (:,:) = ( zua_tlin(:,:,2) + zub_tlin(:,:,2) ) / 10. !-------------------------------------------------------------------- ! Call the tangent routine: dy = L dx !-------------------------------------------------------------------- ua_tl(:,:,:) = zua_tlin(:,:,:) va_tl(:,:,:) = zva_tlin(:,:,:) ub_tl(:,:,:) = zub_tlin(:,:,:) vb_tl(:,:,:) = zvb_tlin(:,:,:) emp_tl (:,:) = zemp_tlin (:,:) sshn_tl(:,:) = zsshn_tlin(:,:) # if defined key_obc gcb_tl (:,:) = zgcb_tlin (:,:) # else gcx_tl (:,:) = 0.e0 ; gcxb_tl(:,:) = 0.e0 gcb_tl (:,:) = 0.e0 gcx_tl (:,:) = zgcx_tlin (:,:) ; gcxb_tl(:,:) = zgcxb_tlin(:,:) # endif CALL dyn_spg_flt_tan( istp, indic ) zua_tlout(:,:,:) = ua_tl(:,:,:) ; zva_tlout(:,:,:) = va_tl(:,:,:) # if defined key_obc zgcx_tlout (:,:) = gcx_tl (:,:) ; zgcxb_tlout(:,:) = gcxb_tl(:,:) # endif zspgu_tlout(:,:) = spgu_tl(:,:) ; zspgv_tlout(:,:) = spgv_tl(:,:) zgcb_tlout (:,:) = gcb_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 # if defined key_obc DO jj = nldj, nlej DO ji = nldi, nlei zgcx_adin (ji,jj) = zgcx_tlout (ji,jj) & & * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,1) * tmask(ji,jj,1) zgcxb_adin(ji,jj) = zgcxb_tlout(ji,jj) & & * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,1) * tmask(ji,jj,1) END DO END DO # endif DO jj = nldj, nlej DO ji = nldi, nlei zgcb_adin (ji,jj) = zgcb_tlout (ji,jj) & & * e1t(ji,jj) * e2t(ji,jj) * fse3u(ji,jj,1) * tmask(ji,jj,1) zspgu_adin (ji,jj) = zspgu_tlout (ji,jj) & & * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1) * umask(ji,jj,1) zspgv_adin(ji,jj) = zspgv_tlout(ji,jj) & & * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1) * vmask(ji,jj,1) END DO END DO !-------------------------------------------------------------------- ! Compute the scalar product: ( L dx )^T W dy !-------------------------------------------------------------------- zsp1 = DOT_PRODUCT( zua_tlout , zua_adin ) & # if defined key_obc & + DOT_PRODUCT( zgcx_tlout , zgcx_adin ) & & + DOT_PRODUCT( zgcxb_tlout, zgcxb_adin ) & # endif & + DOT_PRODUCT( zgcb_tlout , zgcb_adin ) & & + DOT_PRODUCT( zspgu_tlout , zspgu_adin ) & & + DOT_PRODUCT( zspgv_tlout , zspgv_adin ) & & + DOT_PRODUCT( zva_tlout , zva_adin ) !-------------------------------------------------------------------- ! Call the adjoint routine: dx^* = L^T dy^* !-------------------------------------------------------------------- ua_ad(:,:,:) = zua_adin(:,:,:) va_ad(:,:,:) = zva_adin(:,:,:) # if defined key_obc gcx_ad (:,:) = zgcx_adin (:,:) ; gcxb_ad(:,:) = zgcxb_adin(:,:) # else gcx_ad (:,:) = 0.0_wp ; gcxb_ad(:,:) = 0.0_wp gcb_ad (:,:) = zgcb_adin (:,:) spgu_ad(:,:) = zspgu_adin(:,:) spgv_ad(:,:) = zspgv_adin(:,:) ub_ad (:,:,:) = zub_adin (:,:,:) ; vb_ad (:,:,:) = zvb_adin (:,:,:) # endif CALL dyn_spg_flt_adj( istp, indic ) zua_adout(:,:,:) = ua_ad(:,:,:) zva_adout(:,:,:) = va_ad(:,:,:) zub_adout(:,:,:) = ub_ad(:,:,:) zvb_adout(:,:,:) = vb_ad(:,:,:) zemp_adout (:,:) = emp_ad (:,:) zsshn_adout(:,:) = sshn_ad(:,:) zgcx_adout (:,:) = gcx_ad (:,:) zgcxb_adout(:,:) = gcxb_ad(:,:) # if defined key_obc zgcb_adout (:,:) = gcb_ad (:,:) # endif !-------------------------------------------------------------------- ! Compute the scalar product: dx^T L^T W dy !-------------------------------------------------------------------- zsp2 = DOT_PRODUCT( zua_tlin , zua_adout ) & & + DOT_PRODUCT( zva_tlin , zva_adout ) & & + DOT_PRODUCT( zub_tlin , zub_adout ) & & + DOT_PRODUCT( zvb_tlin , zvb_adout ) & & + DOT_PRODUCT( zgcx_tlin , zgcx_adout ) & & + DOT_PRODUCT( zgcxb_tlin, zgcxb_adout ) & # if defined key_obc & + DOT_PRODUCT( zgcb_tlin , zgcb_adout ) & # endif & + DOT_PRODUCT( zsshn_tlin, zsshn_adout ) & & + DOT_PRODUCT( zemp_tlin , zemp_adout ) ! Compare the scalar products ! 14 char:'12345678901234' IF ( istp == nit000 ) THEN SELECT CASE (jpert) CASE(1) cl_name = 'spg_flt Ua T1' CASE(2) cl_name = 'spg_flt Va T1' CASE(3) cl_name = 'spg_flt Ub T1' CASE(4) cl_name = 'spg_flt Vb T1' CASE(5) cl_name = 'spg_flt emp T1' CASE(6) cl_name = 'spg_flt ssh T1' CASE(jpertmax) cl_name = 'dyn_spg_flt T1' END SELECT ELSEIF ( istp == nit000 + 1 ) THEN SELECT CASE (jpert) CASE(1) cl_name = 'spg_flt Ua T2' CASE(2) cl_name = 'spg_flt Va T2' CASE(3) cl_name = 'spg_flt Ub T2' CASE(4) cl_name = 'spg_flt Vb T2' CASE(5) cl_name = 'spg_flt emp T2' CASE(6) cl_name = 'spg_flt ssh T2' CASE(jpertmax) cl_name = 'dyn_spg_flt T2' END SELECT END IF CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) END DO END DO !nn_nmod = kmod ! restore initial frequency of test for the SOR solver ! Deallocate memory DEALLOCATE( & & zua_tlin, & & zva_tlin, & & zub_tlin, & & zvb_tlin, & & zua_tlout, & & zva_tlout, & & zua_adin, & & zva_adin, & & zua_adout, & & zva_adout, & & zub_adout, & & zvb_adout, & & znu & & ) DEALLOCATE( & & zsshn_tlin, & & zemp_tlin, & & zsshn_adout,& & zemp_adout, & & znssh & & ) DEALLOCATE( zgcx_tlin , zgcx_tlout , zgcx_adin , zgcx_adout, & & zgcxb_tlin, zgcxb_tlout, zgcxb_adin, zgcxb_adout, & & zgcb_tlin , zgcb_tlout , zgcb_adin , zgcb_adout & & ) DEALLOCATE ( zub_tlout, zvb_tlout, zub_adin , zvb_adin ) # if defined key_obc DEALLOCATE( zemp_tlout, zemp_adin , & zwn_tlout , zwn_adin ) # endif END SUBROUTINE dyn_spg_flt_adj_tst # if defined key_tst_tlm SUBROUTINE dyn_spg_flt_tlm_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE dyn_spg_flt_tlm_tst *** !! !! ** Purpose : Test the tangent routine. !! !! ** Method : Verify the tangent with Taylor expansion !! !! M(x+hdx) = M(x) + L(hdx) + O(h^2) !! !! where L = tangent routine !! M = direct routine !! dx = input perturbation (random field) !! h = ration on perturbation !! !! In the tangent test we verify that: !! M(x+h*dx) - M(x) !! g(h) = ------------------ ---> 1 as h ---> 0 !! L(h*dx) !! and !! g(h) - 1 !! f(h) = ---------- ---> k (costant) as h ---> 0 !! p !! !! History : !! ! 09-08 (A. Vigilant) !!----------------------------------------------------------------------- !! * Modules used USE dynspg_flt USE tamtrj ! writing out state trajectory USE par_tlm, ONLY: & & tlm_bch, & & cur_loop, & & h_ratio USE istate_mod USE gridrandom, ONLY: & & grid_rd_sd USE trj_tam USE oce , ONLY: & ! ocean dynamics and tracers variables & ua, va, ub, vb, & & sshn USE sbc_oce , ONLY: & & emp USE tamctl, ONLY: & ! Control parameters & numtan, numtan_sc !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit !! * Local declarations REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zua_tlin, & ! Tangent input: ua_tl & zva_tlin, & ! Tangent input: va_tl & zub_tlin, & ! Tangent input: ub_tl & zvb_tlin, & ! & zua_out, & ! & zva_out, & ! & zua_wop, & ! & zva_wop, & & z3r ! 3D field REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & zsshn_tlin, & ! Tangent input: sshn_tl & zemp_tlin, & ! Tangent input: emp_tl & zsshn_out, & ! & zsshn_wop, & & z2r ! 2D field REAL(wp) :: & & zsp1, zsp1_1, zsp1_2, zsp1_3, zsp1_4, & ! & zsp2, zsp2_1, zsp2_2, zsp2_3, zsp2_4, & ! & zsp3, zsp3_1, zsp3_2, zsp3_3, zsp3_4, & ! & zzsp, zzsp_1, zzsp_2, zzsp_3, zzsp_4, & & gamma, & & zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, & & zgsp6, zgsp7 INTEGER :: & & indic, & & istp INTEGER :: & & ji, & & jj, & & jk CHARACTER(LEN=14) :: cl_name CHARACTER (LEN=128) :: file_out, file_wop, file_xdx CHARACTER (LEN=90) :: FMT REAL(KIND=wp), DIMENSION(100):: & & zscua, zscva, & & zscerrua, zscerrva, & & zscsshb, zscsshn, & & zscerrsshb, zscerrsshn INTEGER, DIMENSION(100):: & & iipossshb, iipossshn, iiposua, iiposva, & & ijpossshb, ijpossshn, ijposua, ijposva, & & ikposua, ikposva INTEGER:: & & ii, & & isamp=40, & & jsamp=40, & & ksamp=40, & & numsctlm REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: & & zerrua, zerrva REAL(KIND=wp), DIMENSION(jpi,jpj) :: & & zerrsshb, zerrsshn ! Allocate memory ALLOCATE( & & zua_tlin(jpi,jpj,jpk), & & zva_tlin(jpi,jpj,jpk), & & zub_tlin(jpi,jpj,jpk), & & zvb_tlin(jpi,jpj,jpk), & & zua_out( jpi,jpj,jpk), & & zva_out( jpi,jpj,jpk), & & zua_wop( jpi,jpj,jpk), & & zva_wop( jpi,jpj,jpk), & & z3r(jpi,jpj,jpk) & & ) ALLOCATE( & & zsshn_tlin(jpi,jpj), & & zemp_tlin(jpi,jpj), & & zsshn_out(jpi,jpj), & & zsshn_wop(jpi,jpj), & & z2r(jpi,jpj) & & ) !-------------------------------------------------------------------- ! Reset variables !-------------------------------------------------------------------- zua_tlin( :,:,:) = 0.0_wp zva_tlin( :,:,:) = 0.0_wp zub_tlin( :,:,:) = 0.0_wp zvb_tlin( :,:,:) = 0.0_wp zsshn_tlin( :,:) = 0.0_wp zemp_tlin( :,:) = 0.0_wp zua_out ( :,:,:) = 0.0_wp zva_out ( :,:,:) = 0.0_wp zsshn_out ( :,:) = 0.0_wp zua_wop ( :,:,:) = 0.0_wp zva_wop ( :,:,:) = 0.0_wp zsshn_wop( :,:) = 0.0_wp zscua(:) = 0.0_wp zscva(:) = 0.0_wp zscerrua(:) = 0.0_wp zscerrva(:) = 0.0_wp zerrua(:,:,:) = 0.0_wp zerrva(:,:,:) = 0.0_wp zscsshn(:) = 0.0_wp zscerrsshn(:) = 0.0_wp zerrsshn(:,:) = 0.0_wp !-------------------------------------------------------------------- ! Output filename Xn=F(X0) !-------------------------------------------------------------------- !!!! CALL tlm_namrd gamma = h_ratio file_wop='trj_wop_dynspg' file_xdx='trj_xdx_dynspg' !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN CALL grid_rd_sd( 456953, z3r, 'U', 0.0_wp, stdu) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zua_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 3434334, z3r, 'V', 0.0_wp, stdv) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zva_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 935678, z3r, 'U', 0.0_wp, stdu) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zub_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 1462578, z3r, 'V', 0.0_wp, stdv) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zvb_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 25674910, z2r, 'T', 0.0_wp, stdemp) DO jj = nldj, nlej DO ji = nldi, nlei zemp_tlin(ji,jj) = z2r(ji,jj) END DO END DO CALL grid_rd_sd( 12672456, z2r,'T', 0.0_wp, stdssh) DO jj = nldj, nlej DO ji = nldi, nlei zsshn_tlin(ji,jj) = z2r(ji,jj) END DO END DO ENDIF !-------------------------------------------------------------------- ! Complete Init for Direct !------------------------------------------------------------------- CALL istate_p ! *** initialize the reference trajectory ! ------------ CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN zua_tlin(:,:,:) = gamma * zua_tlin(:,:,:) ua(:,:,:) = ua(:,:,:) + zua_tlin(:,:,:) zva_tlin(:,:,:) = gamma * zva_tlin(:,:,:) va(:,:,:) = va(:,:,:) + zva_tlin(:,:,:) zub_tlin(:,:,:) = gamma * zub_tlin(:,:,:) ub(:,:,:) = ub(:,:,:) + zub_tlin(:,:,:) zvb_tlin(:,:,:) = gamma * zvb_tlin(:,:,:) vb(:,:,:) = vb(:,:,:) + zvb_tlin(:,:,:) zemp_tlin(:,:) = gamma * zemp_tlin(:,:) emp(:,:) = emp(:,:) + zemp_tlin(:,:) zsshn_tlin(:,:) = gamma * zsshn_tlin(:,:) sshn(:,:) = sshn(:,:) + zsshn_tlin(:,:) ENDIF !-------------------------------------------------------------------- ! Compute the direct model F(X0,t=n) = Xn !-------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL dyn_spg_flt(nit000, indic) IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop) IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx) !-------------------------------------------------------------------- ! Compute the Tangent !-------------------------------------------------------------------- IF ( tlm_bch == 2 ) THEN !-------------------------------------------------------------------- ! Initialize the tangent variables !-------------------------------------------------------------------- CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) ua_tl (:,:,:) = zua_tlin (:,:,:) va_tl (:,:,:) = zva_tlin (:,:,:) ub_tl (:,:,:) = zub_tlin (:,:,:) vb_tl (:,:,:) = zvb_tlin (:,:,:) emp_tl (:,: ) = zemp_tlin (:,: ) sshn_tl(:,: ) = zsshn_tlin(:,: ) CALL dyn_spg_flt_tan(nit000, indic) !-------------------------------------------------------------------- ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) !-------------------------------------------------------------------- zsp2_1 = DOT_PRODUCT( ua_tl, ua_tl ) zsp2_2 = DOT_PRODUCT( va_tl, va_tl ) zsp2_4 = DOT_PRODUCT( sshn_tl, sshn_tl ) zsp2 = zsp2_1 + zsp2_2 + zsp2_4 !-------------------------------------------------------------------- ! Storing data !-------------------------------------------------------------------- CALL trj_rd_spl(file_wop) zua_wop (:,:,:) = ua (:,:,:) zva_wop (:,:,:) = va (:,:,:) zsshn_wop(:,:) = sshn(:,:) CALL trj_rd_spl(file_xdx) zua_out (:,:,:) = ua (:,:,:) zva_out (:,:,:) = va (:,:,:) zsshn_out(:,: ) = sshn (:,: ) !-------------------------------------------------------------------- ! Compute the Linearization Error ! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn) ! and ! Compute the Linearization Error ! En = Nn -TL(gamma.dX0, t0,tn) !-------------------------------------------------------------------- ! Warning: Here we re-use local variables z()_out and z()_wop ii=0 DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi zua_out (ji,jj,jk) = zua_out (ji,jj,jk) - zua_wop (ji,jj,jk) zua_wop (ji,jj,jk) = zua_out (ji,jj,jk) - ua_tl (ji,jj,jk) IF ( ua_tl(ji,jj,jk) .NE. 0.0_wp ) & & zerrua(ji,jj,jk) = zua_out(ji,jj,jk)/ua_tl(ji,jj,jk) IF( (MOD(ji, isamp) .EQ. 0) .AND. & & (MOD(jj, jsamp) .EQ. 0) .AND. & & (MOD(jk, ksamp) .EQ. 0) ) THEN ii = ii+1 iiposua(ii) = ji ijposua(ii) = jj ikposua(ii) = jk IF ( INT(umask(ji,jj,jk)) .NE. 0) THEN zscua (ii) = zua_wop(ji,jj,jk) zscerrua (ii) = ( zerrua(ji,jj,jk) - 1.0_wp ) /gamma ENDIF ENDIF END DO END DO END DO ii=0 DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi zva_out (ji,jj,jk) = zva_out (ji,jj,jk) - zva_wop (ji,jj,jk) zva_wop (ji,jj,jk) = zva_out (ji,jj,jk) - va_tl (ji,jj,jk) IF ( va_tl(ji,jj,jk) .NE. 0.0_wp ) & & zerrva(ji,jj,jk) = zva_out(ji,jj,jk)/va_tl(ji,jj,jk) IF( (MOD(ji, isamp) .EQ. 0) .AND. & & (MOD(jj, jsamp) .EQ. 0) .AND. & & (MOD(jk, ksamp) .EQ. 0) ) THEN ii = ii+1 iiposva(ii) = ji ijposva(ii) = jj ikposva(ii) = jk IF ( INT(vmask(ji,jj,jk)) .NE. 0) THEN zscva (ii) = zua_wop(ji,jj,jk) zscerrva (ii) = ( zerrva(ji,jj,jk) - 1.0_wp ) /gamma ENDIF ENDIF END DO END DO END DO ii=0 DO jj = 1, jpj DO ji = 1, jpi zsshn_out (ji,jj) = zsshn_out (ji,jj) - zsshn_wop (ji,jj) zsshn_wop (ji,jj) = zsshn_out (ji,jj) - sshn_tl (ji,jj) IF ( sshn_tl(ji,jj) .NE. 0.0_wp ) & & zerrsshn(ji,jj) = zsshn_out(ji,jj)/sshn_tl(ji,jj) IF( (MOD(ji, isamp) .EQ. 0) .AND. & & (MOD(jj, jsamp) .EQ. 0) ) THEN ii = ii+1 iipossshn(ii) = ji ijpossshn(ii) = jj IF ( INT(tmask(ji,jj,1)) .NE. 0) THEN zscsshn (ii) = zsshn_wop(ji,jj) zscerrsshn (ii) = ( zerrsshn(ji,jj) - 1.0_wp ) /gamma ENDIF ENDIF END DO END DO zsp1_1 = DOT_PRODUCT( zua_out, zua_out ) zsp1_2 = DOT_PRODUCT( zva_out, zva_out ) zsp1_4 = DOT_PRODUCT( zsshn_out, zsshn_out ) zsp1 = zsp1_1 + zsp1_2 + zsp1_4 zsp3_1 = DOT_PRODUCT( zua_wop, zua_wop ) zsp3_2 = DOT_PRODUCT( zva_wop, zva_wop ) zsp3_4 = DOT_PRODUCT( zsshn_wop, zsshn_wop ) zsp3 = zsp3_1 + zsp3_2 + zsp3_4 !-------------------------------------------------------------------- ! Print the linearization error En - norme 2 !-------------------------------------------------------------------- ! 14 char:'12345678901234' cl_name = 'dynspg_tam:En ' zzsp = SQRT(zsp3) zzsp_1 = SQRT(zsp3_1) zzsp_2 = SQRT(zsp3_2) zzsp_4 = SQRT(zsp3_4) zgsp5 = zzsp CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Compute TLM norm2 !-------------------------------------------------------------------- zzsp = SQRT(zsp2) zzsp_1 = SQRT(zsp2_1) zzsp_2 = SQRT(zsp2_2) zzsp_4 = SQRT(zsp2_4) zgsp4 = zzsp cl_name = 'dynspg_tam:Ln2' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Print the linearization error Nn - norme 2 !-------------------------------------------------------------------- zzsp = SQRT(zsp1) zzsp_1 = SQRT(zsp1_1) zzsp_2 = SQRT(zsp1_2) zzsp_4 = SQRT(zsp1_4) cl_name = 'dynspg:Mhdx-Mx' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) zgsp3 = SQRT( zsp3/zsp2 ) zgsp7 = zgsp3/gamma zgsp1 = zzsp zgsp2 = zgsp1 / zgsp4 zgsp6 = (zgsp2 - 1.0_wp)/gamma FMT = "(A8,2X,I4.4,2X,E6.1,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13)" WRITE(numtan,FMT) 'dspgflt ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7 !-------------------------------------------------------------------- ! Unitary calculus !-------------------------------------------------------------------- FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)" cl_name = 'dspgflt ' IF(lwp) THEN DO ii=1, 100, 1 IF ( zscua(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscua ', & & cur_loop, h_ratio, ii, iiposua(ii), ijposua(ii), zscua(ii) ENDDO DO ii=1, 100, 1 IF ( zscva(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscva ', & & cur_loop, h_ratio, ii, iiposva(ii), ijposva(ii), zscva(ii) ENDDO DO ii=1, 100, 1 IF ( zscsshn(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscsshn ', & & cur_loop, h_ratio, ii, iipossshn(ii), ijpossshn(ii), zscsshn(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrua(ii) .NE. 0.0_wp ) WRITE(numsctlm,FMT) cl_name, 'zscerrua ', & & cur_loop, h_ratio, ii, iiposua(ii), ijposua(ii), zscerrua(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrva(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrva ', & & cur_loop, h_ratio, ii, iiposva(ii), ijposva(ii), zscerrva(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrsshn(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrsshn ', & & cur_loop, h_ratio, ii, iipossshn(ii), ijpossshn(ii), zscerrsshn(ii) ENDDO ! write separator WRITE(numtan_sc,"(A4)") '====' ENDIF ENDIF DEALLOCATE( & & zemp_tlin, zsshn_tlin, & & zua_tlin, zva_tlin, & & zua_out, zva_out, & & zua_wop, zva_wop, & & zsshn_out, zsshn_wop, & & z3r, z2r & & ) END SUBROUTINE dyn_spg_flt_tlm_tst # endif # else !!---------------------------------------------------------------------- !! Default case : Empty module No standart explicit free surface !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_spg_flt_tan( kt, kindic ) ! Empty routine WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt END SUBROUTINE dyn_spg_flt_tan SUBROUTINE dyn_spg_flt_adj( kt, kindic ) ! Empty routine WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt END SUBROUTINE dyn_spg_flt_adj SUBROUTINE dyn_spg_flt_adj_tst( kt ) ! Empty routine WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt END SUBROUTINE dyn_spg_flt_adj_tst SUBROUTINE dyn_spg_flt_tlm_tst( kt ) ! Empty routine WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt END SUBROUTINE dyn_spg_flt_tlm_tst # endif #endif END MODULE dynspg_flt_tam