MODULE dynldf_bilap_tam #ifdef key_tam !!=========================================================================== !! *** MODULE dynldf_bilap_tam *** !! Ocean dynamics: lateral viscosity trend !! Tangent and Adjoint Module !!=========================================================================== !!--------------------------------------------------------------------------- !! dyn_ldf_bilap_tan : update the momentum trend with the lateral diffusion !! using an iso-level bilaplacian operator (tangent) !! dyn_ldf_bilap_adj : update the momentum trend with the lateral diffusion !! using an iso-level bilaplacian operator (adjoint) !!--------------------------------------------------------------------------- !! * Modules used USE par_kind , ONLY: & ! Precision variables & wp USE lbclnk , ONLY: & ! Boundary/halo exchange & lbc_lnk USE lbclnk_tam , ONLY: & ! Boundary/halo exchange (adjoint) & lbc_lnk_adj USE par_oce , ONLY: & ! Ocean space and time domain variables & jpi, & & jpj, & & jpk, & & jpim1, & & jpjm1, & & jpkm1 USE oce_tam , ONLY: & & ua_tl, & & va_tl, & & ua_ad, & & va_ad, & & rotb_tl, & & hdivb_tl, & & rotb_ad, & & hdivb_ad USE ldfdyn_oce , ONLY: & ! ocean dynamics: lateral physics & ahm3, & & ahm4, & & ahm0 USE dom_oce , ONLY: & ! Ocean space and time domain & ln_sco, & & ln_zps, & & fmask, & & e1u, & & e2u, & & e1v, & & e2v, & & e1t, & & e2t, & & e1f, & & e2f, & #if defined key_zco & e3t_0 #else & e3u, & & e3v, & & e3t, & & e3f #endif USE in_out_manager, ONLY: & ! I/O manager & lwp, & & numout, & & nit000, & & nitend IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC dyn_ldf_bilap_tan ! called by dynldf_tam.F90 PUBLIC dyn_ldf_bilap_adj ! called by dynldf_tam.F90 !! * Substitutions # include "domzgr_substitute.h90" # include "ldfdyn_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_ldf_bilap_tan( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_bilap_tan *** !! !! ** Purpose : Compute the before trend of the lateral momentum !! diffusion and add it to the general trend of momentum equation. !! !! ** Method : The before horizontal momentum diffusion trend is a !! bi-harmonic operator (bilaplacian type) which separates the !! divergent and rotational parts of the flow. !! Its horizontal components are computed as follow: !! laplacian: !! zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ] !! zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ] !! third derivative: !! * multiply by the eddy viscosity coef. at u-, v-point, resp. !! zlu = ahmu * zlu !! zlv = ahmv * zlv !! * curl and divergence of the laplacian !! zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] ) !! zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] ) !! bilaplacian: !! diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ] !! diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ] !! If ln_sco=F and ln_zps=F, the vertical scale factors in the !! rotational part of the diffusion are simplified !! Add this before trend to the general trend (ua,va): !! (ua,va) = (ua,va) + (diffu,diffv) !! 'key_trddyn' defined: the two components of the horizontal !! diffusion trend are saved. !! !! ** Action : - Update (ua,va) with the before iso-level biharmonic !! mixing trend. !! !! History : !! ! 90-09 (G. Madec) Original code !! ! 91-11 (G. Madec) !! ! 93-03 (M. Guyon) symetrical conditions (M. Guyon) !! ! 96-01 (G. Madec) statement function for e3 !! ! 97-07 (G. Madec) lbc calls !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !! 9.0 ! 04-08 (C. Talandier) New trends organization !! History of the tangent routine !! 9.0 ! 09-12 (F. Vigilant) tangent of 9.0 !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zuatl, zvatl, zbt, ze2u, ze2v ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: & zcutl, zcvtl ! temporary workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zuftl, zuttl, zlutl, zlvtl ! temporary workspace !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap_tan: iso-level bilaplacien operator' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~ ' ENDIF zuftl(:,:,:) = 0.0_wp zuttl(:,:,:) = 0.0_wp zlutl(:,:,:) = 0.0_wp zlvtl(:,:,:) = 0.0_wp ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Laplacian ! --------- IF( ln_sco .OR. ln_zps ) THEN ! s-coordinate or z-coordinate with partial steps zuftl(:,:,jk) = rotb_tl(:,:,jk) * fse3f(:,:,jk) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zlutl(ji,jj,jk) = - ( zuftl(ji,jj,jk) - zuftl(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) & & + ( hdivb_tl(ji+1,jj,jk) - hdivb_tl(ji,jj,jk) ) / e1u(ji,jj) zlvtl(ji,jj,jk) = + ( zuftl(ji,jj,jk) - zuftl(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) & & + ( hdivb_tl(ji,jj+1,jk) - hdivb_tl(ji,jj,jk) ) / e2v(ji,jj) END DO END DO ELSE ! z-coordinate - full step DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zlutl(ji,jj,jk) = - ( rotb_tl (ji ,jj,jk) - rotb_tl (ji,jj-1,jk) ) / e2u(ji,jj) & & + ( hdivb_tl(ji+1,jj,jk) - hdivb_tl(ji,jj ,jk) ) / e1u(ji,jj) zlvtl(ji,jj,jk) = + ( rotb_tl (ji,jj ,jk) - rotb_tl (ji-1,jj,jk) ) / e1v(ji,jj) & & + ( hdivb_tl(ji,jj+1,jk) - hdivb_tl(ji ,jj,jk) ) / e2v(ji,jj) END DO END DO ENDIF ENDDO ! Boundary conditions on the laplacian (zlu,zlv) CALL lbc_lnk( zlutl, 'U', -1.0_wp ) CALL lbc_lnk( zlvtl, 'V', -1.0_wp ) DO jk = 1, jpkm1 ! Third derivative ! ---------------- ! Multiply by the eddy viscosity coef. (at u- and v-points) zlutl(:,:,jk) = zlutl(:,:,jk) * fsahmu(:,:,jk) zlvtl(:,:,jk) = zlvtl(:,:,jk) * fsahmv(:,:,jk) ! Contravariant "laplacian" zcutl(:,:) = e1u(:,:) * zlutl(:,:,jk) zcvtl(:,:) = e2v(:,:) * zlvtl(:,:,jk) ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zuftl(ji,jj,jk) = fmask(ji,jj,jk) * ( zcvtl(ji+1,jj ) - zcvtl(ji,jj) & & - zcutl(ji ,jj+1) + zcutl(ji,jj) ) & #if defined key_zco & / ( e1f(ji,jj)*e2f(ji,jj) ) #else & * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) ) #endif END DO END DO ! Laplacian Horizontal fluxes DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. #if defined key_zco zlutl(ji,jj,jk) = e2u(ji,jj) * zlutl(ji,jj,jk) zlvtl(ji,jj,jk) = e1v(ji,jj) * zlvtl(ji,jj,jk) #else zlutl(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlutl(ji,jj,jk) zlvtl(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlvtl(ji,jj,jk) #endif END DO END DO ! Laplacian divergence DO jj = 2, jpj DO ji = fs_2, jpi ! vector opt. #if defined key_zco zbt = e1t(ji,jj) * e2t(ji,jj) #else zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) #endif zuttl(ji,jj,jk) = ( zlutl(ji,jj,jk) - zlutl(ji-1,jj ,jk) & & + zlvtl(ji,jj,jk) - zlvtl(ji ,jj-1,jk) ) / zbt END DO END DO END DO ! boundary conditions on the laplacian curl and div (zuf,zut) !!bug gm no need to do this 2 following lbc... CALL lbc_lnk( zuftl, 'F', 1.0_wp ) CALL lbc_lnk( zuttl, 'T', 1.0_wp ) DO jk = 1, jpkm1 ! Bilaplacian ! ----------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #if defined key_zco ze2u = e2u(ji,jj) ze2v = e1v(ji,jj) #else ze2u = e2u(ji,jj) * fse3u(ji,jj,jk) ze2v = e1v(ji,jj) * fse3v(ji,jj,jk) #endif ! horizontal biharmonic diffusive trends zuatl = - ( zuftl(ji ,jj,jk) - zuftl(ji,jj-1,jk) ) / ze2u & & + ( zuttl(ji+1,jj,jk) - zuttl(ji,jj ,jk) ) / e1u(ji,jj) zvatl = + ( zuftl(ji,jj ,jk) - zuftl(ji-1,jj,jk) ) / ze2v & & + ( zuttl(ji,jj+1,jk) - zuttl(ji ,jj,jk) ) / e2v(ji,jj) ! add it to the general momentum trends ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + zuatl va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + zvatl END DO END DO ! ! =============== END DO ! End of slab ! ! =============== END SUBROUTINE dyn_ldf_bilap_tan SUBROUTINE dyn_ldf_bilap_adj( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_bilap_adj *** !! !! ** Purpose : Compute the before trend of the lateral momentum !! diffusion and add it to the general trend of momentum equation. !! !! ** Method : The before horizontal momentum diffusion trend is a !! bi-harmonic operator (bilaplacian type) which separates the !! divergent and rotational parts of the flow. !! Its horizontal components are computed as follow: !! laplacian: !! zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ] !! zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ] !! third derivative: !! * multiply by the eddy viscosity coef. at u-, v-point, resp. !! zlu = ahmu * zlu !! zlv = ahmv * zlv !! * curl and divergence of the laplacian !! zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] ) !! zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] ) !! bilaplacian: !! diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ] !! diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ] !! If ln_sco=F and ln_zps=F, the vertical scale factors in the !! rotational part of the diffusion are simplified !! Add this before trend to the general trend (ua,va): !! (ua,va) = (ua,va) + (diffu,diffv) !! 'key_trddyn' defined: the two components of the horizontal !! diffusion trend are saved. !! !! ** Action : - Update (ua,va) with the before iso-level biharmonic !! mixing trend. !! !! History : !! ! 90-09 (G. Madec) Original code !! ! 91-11 (G. Madec) !! ! 93-03 (M. Guyon) symetrical conditions (M. Guyon) !! ! 96-01 (G. Madec) statement function for e3 !! ! 97-07 (G. Madec) lbc calls !! 8.5 ! 02-08 (G. Madec) F90: Free form and module !! 9.0 ! 04-08 (C. Talandier) New trends organization !! History of the adjoint routine !! 9.0 ! 09-12 (F. Vigilant) adjoint of 9.0 !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zuaad, zvaad, zbt, ze2u, ze2v ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: & zcuad, zcvad ! temporary workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zufad, zutad, zluad, zlvad ! temporary workspace !!---------------------------------------------------------------------- IF( kt == nitend ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap_adj: bilaplacien operator' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~ ' ENDIF zuaad = 0.0_wp zvaad = 0.0_wp zufad(:,:,:) = 0.0_wp zutad(:,:,:) = 0.0_wp zluad(:,:,:) = 0.0_wp zlvad(:,:,:) = 0.0_wp zcvad(:,:) = 0.0_wp zcuad(:,:) = 0.0_wp DO jk = 1, jpkm1 ! Bilaplacian ! ----------- DO jj = jpjm1, 2, -1 DO ji = fs_jpim1, fs_2, -1 ! vector opt. #if defined key_zco ze2u = e2u(ji,jj) ze2v = e1v(ji,jj) #else ze2u = e2u(ji,jj) * fse3u(ji,jj,jk) ze2v = e1v(ji,jj) * fse3v(ji,jj,jk) #endif ! add it to the general momentum trends zvaad = zvaad + va_ad(ji,jj,jk) zuaad = zuaad + ua_ad(ji,jj,jk) ! horizontal biharmonic diffusive trends zufad(ji ,jj ,jk) = zufad(ji ,jj ,jk) + zvaad / ze2v zufad(ji-1,jj ,jk) = zufad(ji-1,jj ,jk) - zvaad / ze2v zutad(ji ,jj ,jk) = zutad(ji ,jj ,jk) - zvaad / e2v(ji,jj) zutad(ji ,jj+1,jk) = zutad(ji ,jj+1,jk) + zvaad / e2v(ji,jj) zufad(ji ,jj ,jk) = zufad(ji ,jj ,jk) - zuaad / ze2u zufad(ji ,jj-1,jk) = zufad(ji ,jj-1,jk) + zuaad / ze2u zutad(ji ,jj ,jk) = zutad(ji ,jj ,jk) - zuaad / e1u(ji,jj) zutad(ji+1,jj ,jk) = zutad(ji+1,jj ,jk) + zuaad / e1u(ji,jj) zuaad = 0.0_wp zvaad = 0.0_wp END DO END DO ! ! =============== END DO ! End of slab ! ! =============== ! boundary conditions on the laplacian curl and div (zuf,zut) !!bug gm no need to do this 2 following lbc... CALL lbc_lnk_adj( zutad, 'T', 1.0_wp ) CALL lbc_lnk_adj( zufad, 'F', 1.0_wp ) DO jk = 1, jpkm1 ! Third derivative ! ---------------- ! Laplacian divergence DO jj = jpj, 2, -1 DO ji = jpi, fs_2, -1 ! vector opt. #if defined key_zco zbt = e1t(ji,jj) * e2t(ji,jj) #else zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) #endif zluad(ji ,jj ,jk) = zluad(ji ,jj ,jk) + zutad(ji,jj,jk) / zbt zluad(ji-1,jj ,jk) = zluad(ji-1,jj ,jk) - zutad(ji,jj,jk) / zbt zlvad(ji ,jj ,jk) = zlvad(ji ,jj ,jk) + zutad(ji,jj,jk) / zbt zlvad(ji ,jj-1,jk) = zlvad(ji ,jj-1,jk) - zutad(ji,jj,jk) / zbt zutad(ji,jj,jk) = 0.0_wp END DO END DO ! Laplacian Horizontal fluxes DO jj = jpjm1, 1, -1 DO ji = fs_jpim1, 1, -1 ! vector opt. #if defined key_zco zluad(ji,jj,jk) = e2u(ji,jj) * zluad(ji,jj,jk) zlvad(ji,jj,jk) = e1v(ji,jj) * zlvad(ji,jj,jk) #else zluad(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zluad(ji,jj,jk) zlvad(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlvad(ji,jj,jk) #endif END DO END DO ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps) DO jj = jpjm1, 1, -1 DO ji = fs_jpim1, 1, -1 ! vector opt. #if defined key_zco zufad(ji,jj,jk) = fmask(ji,jj,jk) * zufad(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) ) #else zufad(ji,jj,jk) = fmask(ji,jj,jk) * zufad(ji,jj,jk) * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) ) #endif zcvad(ji ,jj ) = zcvad(ji ,jj ) - zufad(ji,jj,jk) zcvad(ji+1,jj ) = zcvad(ji+1,jj ) + zufad(ji,jj,jk) zcuad(ji ,jj ) = zcuad(ji ,jj ) + zufad(ji,jj,jk) zcuad(ji ,jj+1) = zcuad(ji ,jj+1) - zufad(ji,jj,jk) zufad(ji,jj,jk) = 0.0_wp END DO END DO ! Contravariant "laplacian" DO jj = 1, jpj DO ji = 1, jpi zlvad(ji,jj,jk) = zlvad(ji,jj,jk) + e2v(ji,jj) * zcvad(ji,jj) zluad(ji,jj,jk) = zluad(ji,jj,jk) + e1u(ji,jj) * zcuad(ji,jj) zcvad(ji,jj) = 0.0_wp zcuad(ji,jj) = 0.0_wp END DO END DO ! Multiply by the eddy viscosity coef. (at u- and v-points) zluad(:,:,jk) = zluad(:,:,jk) * fsahmu(:,:,jk) zlvad(:,:,jk) = zlvad(:,:,jk) * fsahmv(:,:,jk) END DO ! Boundary conditions on the laplacian (zlu,zlv) CALL lbc_lnk_adj( zlvad, 'V', -1.0_wp ) CALL lbc_lnk_adj( zluad, 'U', -1.0_wp ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Laplacian ! --------- IF( ln_sco .OR. ln_zps ) THEN ! s-coordinate or z-coordinate with partial steps DO jj = jpjm1, 2, -1 DO ji = fs_jpim1, fs_2, -1 ! vector opt. rotb_ad (ji ,jj ,jk) = rotb_ad (ji ,jj ,jk) & & + zlvad(ji,jj,jk) * fse3f(ji,jj ,jk) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) rotb_ad (ji-1,jj ,jk) = rotb_ad (ji-1,jj ,jk) & & - zlvad(ji,jj,jk) * fse3f(ji-1,jj,jk) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) hdivb_ad(ji ,jj ,jk) = hdivb_ad(ji ,jj ,jk) - zlvad(ji,jj,jk) / e2v(ji,jj) hdivb_ad(ji ,jj+1,jk) = hdivb_ad(ji ,jj+1,jk) + zlvad(ji,jj,jk) / e2v(ji,jj) zlvad(ji,jj,jk) = 0.0_wp rotb_ad (ji ,jj ,jk) = rotb_ad (ji ,jj ,jk) & & - zluad(ji,jj,jk) * fse3f(ji,jj ,jk) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) rotb_ad (ji ,jj-1,jk) = rotb_ad (ji ,jj-1,jk) & & + zluad(ji,jj,jk) * fse3f(ji,jj-1,jk) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) hdivb_ad(ji ,jj ,jk) = hdivb_ad(ji ,jj ,jk) - zluad(ji,jj,jk) / e1u(ji,jj) hdivb_ad(ji+1,jj ,jk) = hdivb_ad(ji+1,jj ,jk) + zluad(ji,jj,jk) / e1u(ji,jj) zluad(ji,jj,jk) = 0.0_wp END DO END DO ! DO jj = 1, jpj ! DO ji = 1, jpi ! rotb_ad(ji,jj,jk) = rotb_ad(ji,jj,jk) + zufad(ji,jj,jk) * fse3f(ji,jj,jk) ! zufad(ji,jj,jk) = 0.0_wp ! END DO ! END DO ELSE ! z-coordinate - full step DO jj = jpjm1, 2, -1 DO ji = fs_jpim1, fs_2, -1 ! vector opt. rotb_ad (ji ,jj ,jk) = rotb_ad (ji ,jj ,jk) + zlvad(ji,jj,jk) / e1v(ji,jj) rotb_ad (ji-1,jj ,jk) = rotb_ad (ji-1,jj ,jk) - zlvad(ji,jj,jk) / e1v(ji,jj) hdivb_ad(ji ,jj ,jk) = hdivb_ad(ji ,jj ,jk) - zlvad(ji,jj,jk) / e2v(ji,jj) hdivb_ad(ji ,jj+1,jk) = hdivb_ad(ji ,jj+1,jk) + zlvad(ji,jj,jk) / e2v(ji,jj) ! zlvad(ji,jj,jk) = 0.0_wp rotb_ad (ji ,jj ,jk) = rotb_ad (ji ,jj ,jk) - zluad(ji,jj,jk) / e2u(ji,jj) rotb_ad (ji ,jj-1,jk) = rotb_ad (ji ,jj-1,jk) + zluad(ji,jj,jk) / e2u(ji,jj) hdivb_ad(ji ,jj ,jk) = hdivb_ad(ji ,jj ,jk) - zluad(ji,jj,jk) / e1u(ji,jj) hdivb_ad(ji+1,jj ,jk) = hdivb_ad(ji+1,jj ,jk) + zluad(ji,jj,jk) / e1u(ji,jj) ! zlvad(ji,jj,jk) = 0.0_wp END DO END DO ENDIF ENDDO END SUBROUTINE dyn_ldf_bilap_adj !!====================================================================== #endif END MODULE dynldf_bilap_tam