MODULE wzvmod_tam #if defined key_tam !!========================================================================== !! *** MODULE wzvmod_tam : TANGENT/ADJOINT OF MODULE wzvmod *** !! !! Ocean diagnostic variable : vertical velocity !! !!========================================================================== !! History of the direct module: !! 5.0 ! 90-10 (C. Levy, G. Madec) Original code !! 7.0 ! 96-01 (G. Madec) Statement function for e3 !! 8.5 ! 02-07 (G. Madec) Free form, F90 !! ! 07-07 (D. Storkey) Zero zhdiv at open boundary (BDY) !! History of the TAM module: !! 7.0 ! 95-01 (F. Van den Berghe) !! 8.0 ! 96-01 (A. Weaver) !! 8.1 ! 98-03 (A. Weaver) !! 8.2 ! 00-08 (A. Weaver) !! 9.0 ! 08-06 (A. Vidard) Skeleton !! 9.0 ! 08-07 (A. Weaver) Tam of the 02-07 version !! 9.0 ! 08-07 (A. Vidard) Tam of the 07-07 version !!---------------------------------------------------------------------- !! wzv_tan : Compute the vertical velocity: tangent routine !! wzv_adj : Compute the vertical velocity: adjoint routine !! wzv_adj_tst : Test of the adjoint routine !!---------------------------------------------------------------------- !! * Modules used 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 in_out_manager, ONLY: & ! I/O manager & lwp, & & numout, & & nit000, & & ln_ctl USE dom_oce , ONLY: & ! Ocean space and time domain & e2u, & & e1v, & & e1t, & & e2t, & # if defined key_vvl & e3t_1, & # else # if defined key_zco & e3t_0, & # else & e3t, & # endif # endif # if defined key_zco ! & e3u_0, & scale factor is identical to e3t_0 ! & e3v_0, & # else & e3u, & & e3v, & # endif & lk_vvl, & & rdt, & & neuler, & & tmask, & & mig, & & mjg, & & nldi, & & nldj, & & nlei, & & nlej USE prtctl , ONLY: & ! Print control & prt_ctl USE domvvl , ONLY: & ! Variable volume & mut USE phycst , ONLY: & ! Physical constants & rauw # if defined key_obc USE obc_par , ONLY: & ! Open boundary conditions & lp_obc_east, & & lp_obc_west, & & lp_obc_north, & & lp_obc_south USE obc_oce , ONLY: & ! Open boundary conditions & nie0p1, & & nie1p1, & & njn0p1, & & njn1p1, & & nje0, & & nje1, & & niw0, & & niw1, & & njw0, & & njw1, & & nin0, & & nin1, & & nis0, & & nis1, & & njs0, & & njs1 # endif USE lbclnk , ONLY: & ! Lateral boundary conditions & lbc_lnk USE lbclnk_tam , ONLY: & ! TAM lateral boundary conditions & lbc_lnk_adj USE divcur_tam , ONLY: & ! TAM horizontal divergence and relative & div_cur_tan ! vorticity USE oce_tam , ONLY: & ! TAM ocean dynamics and tracers variables & un_tl, & & un_ad, & & vn_tl, & & vn_ad, & & wn_tl, & & wn_ad, & & hdivn_tl, & & hdivn_ad, & & sshb_tl, & & sshb_ad USE sbc_oce_tam , ONLY: & ! surface variables & emp_tl, & & emp_ad 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, & & stdssh, & & stdu, & & stdv IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC wzv_tan, & !: tangent routine called by steptan.F90 & wzv_adj, & !: adjoint routine called by stepadj.F90 & wzv_adj_tst !: adjoint test routine !! * Substitutions # include "domzgr_substitute.h90" CONTAINS SUBROUTINE wzv_tan( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE wzv_tan : TANGENT OF wzv *** !! !! ** Purpose of direct routine : !! Compute the now vertical velocity after the array swap !! !! ** Method of direct routine : Using the incompressibility hypothesis, !! the vertical velocity is computed by integrating the horizontal !! divergence from the bottom to the surface. !! The boundary conditions are w=0 at the bottom (no flux) and, !! in regid-lid case, w=0 at the sea surface. !! !! ** action : wn array : the now vertical velocity !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: & & kt ! ocean time-step index !! * Local declarations INTEGER :: & & jk ! dummy loop indices !! Variable volume INTEGER :: & & ji, & ! dummy loop indices & jj REAL(wp) :: & & z2dt, & ! temporary scalar & zraur REAL(wp), DIMENSION (jpi,jpj) :: & & zssha, & & zun, & & zvn, & & zhdiv !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'wzv_tan : vertical velocity from continuity eq.' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ! bottom boundary condition: w=0 (set once for all) wn_tl(:,:,jpk) = 0.0_wp ENDIF IF( lk_vvl ) THEN ! Variable volume ! z2dt = 2.0_wp * rdt ! time step: leap-frog IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest zraur = 1.0_wp / rauw ! Vertically integrated quantities ! -------------------------------- zun(:,:) = 0.0_wp zvn(:,:) = 0.0_wp ! DO jk = 1, jpkm1 ! Vertically integrated transports (now) zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un_tl(:,:,jk) zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn_tl(:,:,jk) END DO ! Horizontal divergence of barotropic transports !-------------------------------------------------- zhdiv(:,:) = 0.0_wp DO jj = 2, jpjm1 DO ji = 2, jpim1 ! vector opt. zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun(ji ,jj ) & & - e2u(ji-1,jj ) * zun(ji-1,jj ) & & + e1v(ji ,jj ) * zvn(ji ,jj ) & & - e1v(ji ,jj-1) * zvn(ji ,jj-1) ) & & / ( e1t(ji,jj) * e2t(ji,jj) ) END DO END DO # if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) ! open boundaries (div must be zero behind the open boundary) ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.0_wp ! east IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.0_wp ! west IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.0_wp ! north IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.0_wp ! south # endif # if defined key_bdy jgrd=1 !: tracer grid. DO jb = 1, nblenrim(jgrd) ji = nbi(jb,jgrd) jj = nbj(jb,jgrd) zhdiv(ji,jj) = 0.e0 END DO # endif CALL lbc_lnk( zhdiv, 'T', 1.0_wp ) ! Sea surface elevation time stepping ! ----------------------------------- zssha(:,:) = sshb_tl(:,:) - z2dt * ( zraur * emp_tl(:,:) & & + zhdiv(:,:) ) * tmask(:,:,1) ! Vertical velocity computed from bottom ! -------------------------------------- DO jk = jpkm1, 1, -1 wn_tl(:,:,jk) = wn_tl(:,:,jk+1) - fse3t(:,:,jk) * hdivn_tl(:,:,jk) & & - ( zssha(:,:) - sshb_tl(:,:) ) & & * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt END DO ELSE ! Fixed volume ! Vertical velocity computed from bottom ! -------------------------------------- DO jk = jpkm1, 1, -1 wn_tl(:,:,jk) = wn_tl(:,:,jk+1) - fse3t(:,:,jk) * hdivn_tl(:,:,jk) END DO ENDIF IF(ln_ctl) CALL prt_ctl(tab3d_1=wn_tl, clinfo1=' w**2 - : ', mask1=wn_tl) END SUBROUTINE wzv_tan SUBROUTINE wzv_adj( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE wzv_adj : ADJOINT OF wzv_tan *** !! !! ** Purpose of direct routine : !! Compute the now vertical velocity after the array swap !! !! ** Method of direct routine : Using the incompressibility hypothesis, !! the vertical velocity is computed by integrating the horizontal !! divergence from the bottom to the surface. !! The boundary conditions are w=0 at the bottom (no flux) and, !! in regid-lid case, w=0 at the sea surface. !! !! ** action : wn array : the now vertical velocity !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local declarations INTEGER :: & & jk ! dummy loop indices !! Variable volume INTEGER :: & & ji, & ! dummy loop indices & jj REAL(wp) :: & & z2dt, & ! temporary scalars & zraur REAL(wp), DIMENSION (jpi,jpj) :: & & zssha, & ! workspace & zun, & & zvn, & & zhdiv, & & zfac1, & & zfac2 IF( lk_vvl ) THEN ! Variable volume ! z2dt = 2.0_wp * rdt ! time step: leap-frog IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest zraur = 1.0_wp / rauw ! Local adjoint variable initialization ! ------------------------------------- zssha(:,:) = 0.0_wp zun (:,:) = 0.0_wp zvn (:,:) = 0.0_wp ! Vertical velocity computed from bottom ! -------------------------------------- DO jk = 1, jpkm1 zfac1(:,:) = fsve3t(:,:,jk) * mut(:,:,jk) / z2dt hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) - wn_ad(:,:,jk) * fse3t(:,:,jk) wn_ad(:,:,jk+1) = wn_ad(:,:,jk+1) + wn_ad(:,:,jk) wn_ad(:,:,jk ) = wn_ad(:,:,jk ) * zfac1(:,:) sshb_ad(:,:) = sshb_ad(:,:) + wn_ad(:,:,jk) zssha(:,:) = zssha(:,:) - wn_ad(:,:,jk) wn_ad(:,:,jk ) = 0.0_wp END DO ! Sea surface elevation time stepping ! ----------------------------------- zfac2(:,:) = z2dt * tmask(:,:,1) zhdiv(:,:) = - zssha(:,:) * zfac2(:,:) emp_ad(:,:) = emp_ad(:,:) - zssha(:,:) * zraur * zfac2(:,:) sshb_ad(:,:) = sshb_ad(:,:) + zssha(:,:) zssha(:,:) = 0.0_wp CALL lbc_lnk_adj( zhdiv, 'T', 1.0_wp ) # if defined key_bdy jgrd=1 !: tracer grid. DO jb = 1, nblenrim(jgrd) ji = nbi(jb,jgrd) jj = nbj(jb,jgrd) zhdiv(ji,jj) = 0.0_wp END DO # endif # if defined key_obc && ( key_dynspg_exp || key_dynspg_ts ) ! open boundaries (div must be zero behind the open boundary) ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.0_wp ! east IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.0_wp ! west IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.0_wp ! north IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.0_wp ! south # endif ! Horizontal divergence of barotropic transports !-------------------------------------------------- DO jj = jpjm1, 2, -1 DO ji = jpim1, 2, -1 ! vector opt. zhdiv(ji,jj) = zhdiv(ji,jj) / ( e1t(ji,jj) * e2t(ji,jj) ) zun(ji ,jj ) = zun(ji ,jj ) + zhdiv(ji,jj) * e2u(ji ,jj ) zun(ji-1,jj ) = zun(ji-1,jj ) - zhdiv(ji,jj) * e2u(ji-1,jj ) zvn(ji ,jj ) = zvn(ji ,jj ) + zhdiv(ji,jj) * e1v(ji ,jj ) zvn(ji ,jj-1) = zvn(ji ,jj-1) - zhdiv(ji,jj) * e1v(ji ,jj-1) zhdiv(ji,jj) = 0.0_wp END DO END DO DO jk = jpkm1, 1, -1 ! Vertically integrated transports (now) un_ad(:,:,jk) = un_ad(:,:,jk) + zun(:,:) * fse3u(:,:,jk) vn_ad(:,:,jk) = vn_ad(:,:,jk) + zvn(:,:) * fse3v(:,:,jk) END DO ! Vertically integrated quantities ! -------------------------------- zun(:,:) = 0.0_wp zvn(:,:) = 0.0_wp ELSE ! Fixed volume ! Vertical velocity computed from bottom ! -------------------------------------- DO jk = 1, jpkm1 hdivn_ad(:,:,jk) = hdivn_ad(:,:,jk) & & - fse3t(:,:,jk) * wn_ad(:,:,jk) wn_ad(:,:,jk+1) = wn_ad(:,:,jk+1) + wn_ad(:,:,jk) wn_ad(:,:,jk ) = 0.0_wp END DO ENDIF END SUBROUTINE wzv_adj SUBROUTINE wzv_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE wzv_adj_tst : TEST OF wzv_adj *** !! !! ** 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: !! !! If variable volume ( lk_vvl = .TRUE ) then !! 1) dx = ( un_tl, vn_tl, emp_tl, sshb_tl ) and !! dy = ( wn_tl ) !! Otherwise !! 2) dx = ( hdivn_tl ) and !! dy = ( wn_tl ) !! !! History : !! ! 08-07 (A. Weaver) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit !! * Local declarations REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zhdivn_tlin, & ! Tangent input: now horizontal divergence & zun_tlin, & ! Tangent input: now u-velocity & zvn_tlin, & ! Tangent input: now v-velocity & zwn_tlout, & ! Tangent output: now w-velocity & zwn_adin, & ! Adjoint input: now w-velocity & zhdivn_adout, & ! Adjoint output: now horizontal divergence & zun_adout, & ! Adjoint output: now u-velocity & zvn_adout, & ! Adjoint output: now v-velocity & znu, & ! 3D random field for u & znv ! 3D random field for v REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & & zsshb_tlin, & ! Tangent input: before SSH & zsshb_adout, & ! Adjoint output: before SSH & zemp_tlin, & ! Tangent input: EmP & zemp_adout, & ! Adjoint output: EmP & znssh, & ! 2D random field for SSH & znemp ! 2D random field for EmP INTEGER :: & & ji, & ! dummy loop indices & jj, & & jk INTEGER, DIMENSION(jpi,jpj) :: & & iseed_2d ! 2D seed for the random number generator REAL(KIND=wp) :: & ! random field standard deviation for: & zstdu, & ! u-velocity & zstdv, & ! v-velocity & zstdssh, & ! SSH & zstdemp, & ! EMP & zsp1, & ! scalar product involving the tangent routine & zsp2, & ! scalar product involving the adjoint routine & zsp2_1, & ! scalar product components & zsp2_2, & & zsp2_3, & & zsp2_4, & & zsp2_5, & & z2dt, & ! temporary scalars & zraur CHARACTER (LEN=14) :: & & cl_name ! Allocate memory ALLOCATE( & & zhdivn_tlin(jpi,jpj,jpk), & & zun_tlin(jpi,jpj,jpk), & & zvn_tlin(jpi,jpj,jpk), & & zwn_tlout(jpi,jpj,jpk), & & zwn_adin(jpi,jpj,jpk), & & zhdivn_adout(jpi,jpj,jpk), & & zun_adout(jpi,jpj,jpk), & & zvn_adout(jpi,jpj,jpk), & & znu(jpi,jpj,jpk), & & znv(jpi,jpj,jpk) & & ) ALLOCATE( & & zsshb_tlin(jpi,jpj), & & zsshb_adout(jpi,jpj), & & zemp_tlin(jpi,jpj), & & zemp_adout(jpi,jpj), & & znssh(jpi,jpj), & & znemp(jpi,jpj) & & ) ! Initialize constants z2dt = 2.0_wp * rdt ! time step: leap-frog zraur = 1.0_wp / rauw ! inverse density of pure water (m3/kg) !============================================================= ! 1) dx = ( un_tl, vn_tl, emp_tl, sshb_tl ) and dy = ( wn_tl ) ! - or - ! 2) dx = ( hdivn_tl ) and dy = ( wn_tl ) !============================================================= !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 785483 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znssh, 'T', 0.0_wp, stdssh ) DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 358606 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, znemp, 'T', 0.0_wp, stdssh ) 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 jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei un_tl(ji,jj,jk) = znu(ji,jj,jk) vn_tl(ji,jj,jk) = znv(ji,jj,jk) END DO END DO END DO CALL div_cur_tan( nit000 ) ! Generate noise hdivn 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) zhdivn_tlin(ji,jj,jk) = hdivn_tl(ji,jj,jk) END DO END DO END DO DO jj = nldj, nlej DO ji = nldi, nlei zsshb_tlin(ji,jj) = znssh(ji,jj) zemp_tlin (ji,jj) = znemp(ji,jj) / ( z2dt * zraur ) END DO END DO !-------------------------------------------------------------------- ! Call the tangent routine: dy = L dx !-------------------------------------------------------------------- un_tl (:,:,:) = zun_tlin (:,:,:) vn_tl (:,:,:) = zvn_tlin (:,:,:) hdivn_tl(:,:,:) = zhdivn_tlin(:,:,:) sshb_tl(:,:) = zsshb_tlin(:,:) emp_tl (:,:) = zemp_tlin (:,:) CALL wzv_tan( nit000 ) zwn_tlout(:,:,:) = wn_tl(:,:,:) !-------------------------------------------------------------------- ! Initialize the adjoint variables: dy^* = W dy !-------------------------------------------------------------------- DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zwn_adin(ji,jj,jk) = zwn_tlout(ji,jj,jk) & & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) & & * tmask(ji,jj,jk) END DO END DO END DO !-------------------------------------------------------------------- ! Compute the scalar product: ( L dx )^T W dy !-------------------------------------------------------------------- zsp1 = DOT_PRODUCT( zwn_tlout, zwn_adin ) !-------------------------------------------------------------------- ! Call the adjoint routine: dx^* = L^T dy^* !-------------------------------------------------------------------- wn_ad(:,:,:) = zwn_adin(:,:,:) CALL wzv_adj( nit000 ) zun_adout (:,:,:) = un_ad (:,:,:) zvn_adout (:,:,:) = vn_ad (:,:,:) zhdivn_adout(:,:,:) = hdivn_ad(:,:,:) zsshb_adout(:,:) = sshb_ad(:,:) zemp_adout (:,:) = emp_ad (:,:) !-------------------------------------------------------------------- ! Compute the scalar product: dx^T L^T W dy !-------------------------------------------------------------------- zsp2_1 = DOT_PRODUCT( zun_tlin, zun_adout ) zsp2_2 = DOT_PRODUCT( zvn_tlin, zvn_adout ) zsp2_3 = DOT_PRODUCT( zhdivn_tlin, zhdivn_adout ) zsp2_4 = DOT_PRODUCT( zemp_tlin, zemp_adout ) zsp2_5 = DOT_PRODUCT( zsshb_tlin, zsshb_adout ) IF( lk_vvl ) THEN zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 ELSE zsp2 = zsp2_3 ENDIF ! Compare the scalar products ! 14 char:'12345678901234' cl_name = 'wzv_adj ' CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) END SUBROUTINE wzv_adj_tst !!====================================================================== #endif END MODULE wzvmod_tam