MODULE eosbn2_tam !!============================================================================== !! *** MODULE eosbn2_tam *** !! Ocean diagnostic variable : equation of state - in situ and potential density !! - Brunt-Vaisala frequency !! Tangent and Adjoint Module !!=========================================================================== !! History of the direct Module : !! OPA ! 1989-03 (O. Marti) Original code !! 6.0 ! 1994-07 (G. Madec, M. Imbard) add bn2 !! 6.0 ! 1994-08 (G. Madec) Add Jackett & McDougall eos !! 7.0 ! 1996-01 (G. Madec) statement function for e3 !! 8.1 ! 1997-07 (G. Madec) density instead of volumic mass !! - ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure gradient !! 8.2 ! 2001-09 (M. Ben Jelloul) bugfix on linear eos !! NEMO 1.0 ! 2002-10 (G. Madec) add eos_init !! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d !! - ! 2003-08 (G. Madec) F90, free form !! 3.0 ! 2006-08 (G. Madec) add tfreez function !! ! 2003-09 (M. Balmaseda) compute refrence rho prof !! History of the TAM Module : !! 8.2 ! 2005-03 (F. Van den Berghe, A. Weaver, N. Daget) !! - eostan.F !! 9.0 ! 2007-07 (K. Mogensen) Initial version based on eostan.F !! ! 2008-07 (A. Vidard) bug fix in computation of prd_tl if neos=1 !! ! 2008-11 (A. Vidard) TAM of the 06-08 version !! NEMO 3.2 ! 2010-04 (F. Vigilant) version 3.2 !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! Direct subroutines !! eos : generic interface of the equation of state !! eos_insitu : Compute the in situ density !! eos_insitu_pot : Compute the insitu and surface referenced potential !! volumic mass !! eos_insitu_2d : Compute the in situ density for 2d fields !! eos_insitu_pot_1pt : Compute the in situ density for a single point !! eos_bn2 : Compute the Brunt-Vaisala frequency !! tfreez : Compute the surface freezing temperature (NOT IN TAM) !! eos_init : set eos parameters (namelist) !! eos_rprof : Compute the in situ density of a reference profile !!---------------------------------------------------------------------- !! * Modules used #if defined key_zdfddm USE oce_tam , ONLY: & & rrau_tl, & & rrau_ad #endif USE dom_oce , ONLY: & ! ocean space and time domain & tmask, & & e1t, & & e2t, & #if defined key_zco & e3t_0, & & e3w_0, & & gdept_0, & & gdepw_0, & #else & e3t, & & e3w, & & gdept, & & gdepw, & #endif & mig, & & mjg, & & nperio, & & nldi, & & nldj, & & nlei, & & nlej USE par_kind , ONLY: & & wp USE par_oce , ONLY: & & jpi, & & jpj, & & jpk, & & jpim1, & & jpjm1, & & jpkm1, & & jpiglo USE oce , ONLY: & & tn, & & sn USE phycst , ONLY: & ! physical constants & rau0, & & grav USE in_out_manager, ONLY: & ! I/O manager & lwp, & & ctmp1, & & numout, & & ctl_stop #if defined key_zdfddm USE zdfddm ! vertical physics: double diffusion #endif USE eosbn2 , ONLY: & & eos_init, & & neos_init, & & nn_eos, & & rn_alpha, & & rn_beta, & & ralpbet 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, & & prntst_tlm, & ! & stds, & & stdt IMPLICIT NONE PRIVATE !! * Interface INTERFACE eos_tan MODULE PROCEDURE eos_insitu_tan, eos_insitu_pot_tan, eos_pot_1pt_tan, eos_insitu_2d_tan END INTERFACE INTERFACE eos_adj MODULE PROCEDURE eos_insitu_adj, eos_insitu_pot_adj, eos_pot_1pt_adj, eos_insitu_2d_adj END INTERFACE INTERFACE bn2_tan MODULE PROCEDURE eos_bn2_tan END INTERFACE INTERFACE bn2_adj MODULE PROCEDURE eos_bn2_adj END INTERFACE !! * Routine accessibility PUBLIC eos_tan ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90 PUBLIC bn2_tan ! called by step.F90 PUBLIC eos_adj ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90 PUBLIC bn2_adj ! called by step.F90 #if defined key_tam PUBLIC eos_adj_tst PUBLIC bn2_adj_tst #if defined key_tst_tlm PUBLIC eos_tlm_tst PUBLIC bn2_tlm_tst #endif #endif !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" CONTAINS SUBROUTINE eos_insitu_tan( ptem, psal, ptem_tl, psal_tl, prd_tl ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_insitu_tan : TL OF ROUTINE eos_insitu *** !! !! ** Purpose of direct routine : Compute the in situ density !! (ratio rho/rau0) from potential temperature and salinity !! using an equation of state defined through the namelist !! parameter nn_eos. !! !! ** Method of direct routine : 3 cases: !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. !! the in situ density is computed directly as a function of !! potential temperature relative to the surface (the opa t !! variable), salt and pressure (assuming no pressure variation !! along geopotential surfaces, i.e. the pressure p in decibars !! is approximated by the depth in meters. !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 !! with pressure p decibars !! potential temperature t deg celsius !! salinity s psu !! reference volumic mass rau0 kg/m**3 !! in situ volumic mass rho kg/m**3 !! in situ density anomalie prd no units !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, !! t = 40 deg celcius, s=40 psu !! nn_eos = 1 : linear equation of state function of temperature only !! prd(t) = 0.0285 - rn_alpha * t !! nn_eos = 2 : linear equation of state function of temperature and !! salinity !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. !! Note that no boundary condition problem occurs in this routine !! as (ptem,psal) are defined over the whole domain. !! !! ** Comments on Adjoint Routine : !! Care has been taken to avoid division by zero when computing !! the inverse of the square root of salinity at masked salinity !! points. !! !! * Arguments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & & ptem, & ! potential temperature & psal, & ! salinity & ptem_tl, & ! TL of potential temperature & psal_tl ! TL of salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: & & prd_tl ! TL of potential density (surface referenced) !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & ! temporary scalars zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, & zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, & zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, & zr4tl, zrhoptl, zetl, zbwtl, & zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, & zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, & zmask, zrau0r REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws !!---------------------------------------------------------------------- SELECT CASE ( nn_eos ) CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==! zrau0r = 1.e0 / rau0 #ifdef key_sp zeps = 1.e-7 #else zeps = 1.e-14 #endif !CDIR NOVERRCHK zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) ) ! DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zt = ptem(ji,jj,jk) zs = psal(ji,jj,jk) zh = fsdept(ji,jj,jk) ! depth zsr= zws(ji,jj,jk) ! square root salinity ! compute volumic mass pure water at atm pressure zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 ! seawater volumic mass atm pressure zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt & -4.0899e-3 ) *zt+0.824493 zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 zr4= 4.8314e-4 ! potential volumic mass (reference to the surface) zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 ! add the compression terms ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 zbw= ( 1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 zb = zbw + ze * zs zd = -2.042967e-2 zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 za = ( zd*zsr + zc ) *zs + zaw zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545 za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 zk0= ( zb1*zsr + za1 )*zs + zkw ! Tangent linear part zttl = ptem_tl(ji,jj,jk) zstl = psal_tl(ji,jj,jk) zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) & & * tmask(ji,jj,jk) * zstl zr1tl= ( ( ( ( 5.*6.536332e-9 * zt & & -4.*1.120083e-6 ) * zt & & +3.*1.001685e-4 ) * zt & & -2.*9.095290e-3 ) * zt & & + 6.793952e-2 ) * zttl zr2tl= ( ( ( 4.*5.3875e-9 * zt & & -3.*8.2467e-7 ) * zt & & +2.*7.6438e-5 ) * zt & & - 4.0899e-3 ) * zttl zr3tl= ( -2.*1.6546e-6 * zt & & + 1.0227e-4 ) * zttl zrhoptl= zr1tl & & + zs * zr2tl & & + zsr * zs * zr3tl & & + zr3 * zs * zsrtl & & + ( 2. * zr4 * zs + zr2 & & + zr3 * zsr ) * zstl zetl = ( -2.*3.508914e-8 * zt & & - 1.248266e-8 ) * zttl zbwtl= ( 2.*1.296821e-6 * zt & & - 5.782165e-9 ) * zttl zbtl = zbwtl & & + zs * zetl & & + ze * zstl zctl = ( -2.*7.267926e-5 * zt & & + 2.598241e-3 ) * zttl zawtl= ( ( 3.*5.939910e-6 * zt & & +2.*2.512549e-3 ) * zt & & - 0.1028859 ) * zttl zatl = zawtl & & + zd * zs * zsrtl & & + zs * zctl & & + ( zd * zsr + zc ) * zstl zb1tl= ( -2.*0.1909078 * zt & & + 7.390729 ) * zttl za1tl= ( ( 3.*2.326469e-3 * zt & & +2.*1.553190 ) * zt & & - 65.00517 ) * zttl zkwtl= ( ( ( -4.*1.361629e-4 * zt & & -3.*1.852732e-2 ) * zt & & -2.*30.41638 ) * zt & & + 2098.925 ) * zttl zk0tl= zkwtl & & + zb1 * zs * zsrtl & & + zs * zsr * zb1tl & & + zs * za1tl & & + ( zb1 * zsr + za1 ) * zstl ! Masked in situ density anomaly zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) ) zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 ) prd_tl(ji,jj,jk) = tmask(ji,jj,jk) * zrdc2 * & & ( zrhoptl & & - zrdc2 * zh * zrdc1**2 * zrhop & & * ( zk0tl & & - zh * ( zatl & & - zh * zbtl ) ) )& & * zrau0r END DO END DO END DO ! CASE ( 1 ) !== Linear formulation function of temperature only ==! DO jk = 1, jpkm1 prd_tl(:,:,jk) = ( - rn_alpha * ptem_tl(:,:,jk) ) * tmask(:,:,jk) END DO ! CASE ( 2 ) !== Linear formulation function of temperature and salinity ==! ! ! =============== DO jk = 1, jpkm1 prd_tl(:,:,jk) = ( rn_beta * psal_tl(:,:,jk) - rn_alpha * ptem_tl(:,:,jk ) ) * tmask(:,:,jk) END DO ! END SELECT END SUBROUTINE eos_insitu_tan SUBROUTINE eos_insitu_pot_tan( ptem, psal, ptem_tl, psal_tl, prd_tl, prhop_tl) !!---------------------------------------------------------------------- !! *** ROUTINE eos_insitu_pot_tan *** !! !! ** Purpose or the direct routine: !! Compute the in situ density (ratio rho/rau0) and the !! potential volumic mass (Kg/m3) from potential temperature and !! salinity fields using an equation of state defined through the !! namelist parameter nn_eos. !! !! ** Method : !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. !! the in situ density is computed directly as a function of !! potential temperature relative to the surface (the opa t !! variable), salt and pressure (assuming no pressure variation !! along geopotential surfaces, i.e. the pressure p in decibars !! is approximated by the depth in meters. !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 !! rhop(t,s) = rho(t,s,0) !! with pressure p decibars !! potential temperature t deg celsius !! salinity s psu !! reference volumic mass rau0 kg/m**3 !! in situ volumic mass rho kg/m**3 !! in situ density anomalie prd no units !! !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, !! t = 40 deg celcius, s=40 psu !! !! nn_eos = 1 : linear equation of state function of temperature only !! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t !! rhop(t,s) = rho(t,s) !! !! nn_eos = 2 : linear equation of state function of temperature and !! salinity !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 !! = rn_beta * s - rn_alpha * tn - 1. !! rhop(t,s) = rho(t,s) !! Note that no boundary condition problem occurs in this routine !! as (tn,sn) or (ta,sa) are defined over the whole domain. !! !! ** Action : - prd , the in situ density (no units) !! - prhop, the potential volumic mass (Kg/m3) !! !! References : !! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994 !! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978 !! !!---------------------------------------------------------------------- !! * Arguments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & ptem, & ! potential temperature psal ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & ptem_tl,& ! potential temperature psal_tl ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: & prd_tl, & ! potential density (surface referenced) prhop_tl ! potential density (surface referenced) !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & ! temporary scalars zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, & zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, & zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, & zr4tl, zrhoptl, zetl, zbwtl, & zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, & zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, & zmask, zrau0r REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws !!---------------------------------------------------------------------- SELECT CASE ( nn_eos ) CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==! zrau0r = 1.e0 / rau0 #ifdef key_sp zeps = 1.e-7 #else zeps = 1.e-14 #endif !CDIR NOVERRCHK zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) ) ! DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zt = ptem(ji,jj,jk) zs = psal(ji,jj,jk) zh = fsdept(ji,jj,jk) ! depth zsr = zws(ji,jj,jk) ! square root salinity ! compute volumic mass pure water at atm pressure zr1 = ( ( ( ( 6.536332e-9 * zt - 1.120083e-6 ) * zt & & + 1.001685e-4 ) * zt - 9.095290e-3 ) * zt & & + 6.793952e-2 ) * zt + 999.842594 ! seawater volumic mass atm pressure zr2 = ( ( ( 5.3875e-9 * zt - 8.2467e-7 ) * zt & & + 7.6438e-5 ) * zt - 4.0899e-3 ) * zt & & + 0.824493 zr3 = ( -1.6546e-6 * zt + 1.0227e-4 ) * zt - 5.72466e-3 zr4 = 4.8314e-4 ! potential volumic mass (reference to the surface) zrhop= ( zr4 * zs + zr3 * zsr + zr2 ) * zs + zr1 ! add the compression terms ze = ( -3.508914e-8 * zt - 1.248266e-8 ) * zt - 2.595994e-6 zbw = ( 1.296821e-6 * zt - 5.782165e-9 ) * zt + 1.045941e-4 zb = zbw + ze * zs zd = -2.042967e-2 zc = (-7.267926e-5 * zt + 2.598241e-3 ) * zt + 0.1571896 zaw= ( ( 5.939910e-6 * zt + 2.512549e-3 ) * zt - 0.1028859 ) * zt - 4.721788 za = ( zd * zsr + zc ) * zs + zaw zb1 = (-0.1909078 * zt + 7.390729 ) * zt - 55.87545 za1 = ( ( 2.326469e-3 * zt + 1.553190 ) * zt - 65.00517 & & ) * zt + 1044.077 zkw = ( ( (-1.361629e-4 * zt - 1.852732e-2 ) * zt - 30.41638 & & ) * zt + 2098.925 ) * zt + 190925.6 zk0 = ( zb1 * zsr + za1 ) * zs + zkw zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) ) zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 ) ! Tangent linear part zttl = ptem_tl(ji,jj,jk) zstl = psal_tl(ji,jj,jk) zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) & & * tmask(ji,jj,jk) * zstl zr1tl= ( ( ( ( 5.*6.536332e-9 * zt & & -4.*1.120083e-6 ) * zt & & +3.*1.001685e-4 ) * zt & & -2.*9.095290e-3 ) * zt & & + 6.793952e-2 ) * zttl zr2tl= ( ( ( 4.*5.3875e-9 * zt & & -3.*8.2467e-7 ) * zt & & +2.*7.6438e-5 ) * zt & & - 4.0899e-3 ) * zttl zr3tl= ( -2.*1.6546e-6 * zt & & + 1.0227e-4 ) * zttl zrhoptl= zr1tl & & + zs * zr2tl & & + zsr * zs * zr3tl & & + zr3 * zs * zsrtl & & + ( 2. * zr4 * zs + zr2 & & + zr3 * zsr ) * zstl prhop_tl(ji,jj,jk) = zrhoptl * tmask(ji,jj,jk) zetl = ( -2.*3.508914e-8 * zt & & - 1.248266e-8 ) * zttl zbwtl= ( 2.*1.296821e-6 * zt & & - 5.782165e-9 ) * zttl zbtl = zbwtl & & + zs * zetl & & + ze * zstl zctl = ( -2.*7.267926e-5 * zt & & + 2.598241e-3 ) * zttl zawtl= ( ( 3.*5.939910e-6 * zt & & +2.*2.512549e-3 ) * zt & & - 0.1028859 ) * zttl zatl = zawtl & & + zd * zs * zsrtl & & + zs * zctl & & + ( zd * zsr + zc ) * zstl zb1tl= ( -2.*0.1909078 * zt & & + 7.390729 ) * zttl za1tl= ( ( 3.*2.326469e-3 * zt & & +2.*1.553190 ) * zt & & - 65.00517 ) * zttl zkwtl= ( ( ( -4.*1.361629e-4 * zt & & -3.*1.852732e-2 ) * zt & & -2.*30.41638 ) * zt & & + 2098.925 ) * zttl zk0tl= zkwtl & & + zb1 * zs * zsrtl & & + zs * zsr * zb1tl & & + zs * za1tl & & + ( zb1 * zsr + za1 ) * zstl ! Masked in situ density anomaly prd_tl(ji,jj,jk) = tmask(ji,jj,jk) * zrdc2 * & & ( zrhoptl & & - zrdc2 * zh * zrdc1**2 * zrhop & & * ( zk0tl & & - zh * ( zatl & & - zh * zbtl ) ) )& & * zrau0r END DO END DO END DO ! CASE ( 1 ) !== Linear formulation = F( temperature ) ==! DO jk = 1, jpkm1 prd_tl (:,:,jk) = ( - rn_alpha * ptem_tl(:,:,jk) ) * tmask(:,:,jk) prhop_tl(:,:,jk) = ( rau0 * prd_tl(:,:,jk) ) * tmask(:,:,jk) END DO ! CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! DO jk = 1, jpkm1 prd_tl(:,:,jk) = ( rn_beta * psal_tl(:,:,jk) - rn_alpha * ptem_tl(:,:,jk) ) * tmask(:,:,jk) prhop_tl(:,:,jk) = ( rau0 * prd_tl(:,:,jk) ) * tmask(:,:,jk) END DO ! END SELECT END SUBROUTINE eos_insitu_pot_tan SUBROUTINE eos_pot_1pt_tan( ptem, psal, ptem_tl, psal_tl, prhop_tl) !!---------------------------------------------------------------------- !! *** ROUTINE eos_pot_1pt_tan *** !! !! ** Purpose of the direct routine: !! Compute the !! potential volumic mass (Kg/m3) from potential temperature and !! salinity fields using an equation of state defined through the !! namelist parameter neos. !! !! ** Method : !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. !! the in situ density is computed directly as a function of !! potential temperature relative to the surface (the opa t !! variable), salt and pressure (assuming no pressure variation !! along geopotential surfaces, i.e. the pressure p in decibars !! is approximated by the depth in meters. !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 !! rhop(t,s) = rho(t,s,0) !! with pressure p decibars !! potential temperature t deg celsius !! salinity s psu !! reference volumic mass rau0 kg/m**3 !! in situ volumic mass rho kg/m**3 !! in situ density anomalie prd no units !! !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, !! t = 40 deg celcius, s=40 psu !! !! nn_eos = 1 : linear equation of state function of temperature only !! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - ralpha * t !! rhop(t,s) = rho(t,s) !! !! nn_eos = 2 : linear equation of state function of temperature and !! salinity !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 !! = rn_beta * s - rn_alpha * tn - 1. !! rhop(t,s) = rho(t,s) !! Note that no boundary condition problem occurs in this routine !! as (tn,sn) or (ta,sa) are defined over the whole domain. !! !! ** Action : - prd , the in situ density (no units) !! - prhop, the potential volumic mass (Kg/m3) !! !! References : !! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994 !! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978 !! !! History of the direct routine: !! 4.0 ! 89-03 (O. Marti) !! ! 94-08 (G. Madec) !! ! 96-01 (G. Madec) statement function for e3 !! ! 97-07 (G. Madec) introduction of neos, OPA8.1 !! ! 97-07 (G. Madec) density instead of volumic mass !! ! 99-02 (G. Madec, N. Grima) semi-implicit pressure gradient !! ! 01-09 (M. Ben Jelloul) bugfix !! 9.0 ! 03-08 (G. Madec) F90, free form !! ! 06-11 (G. Smith) single point case !! History of the tangent routine: !! 9.0 ! 08-07 (A. Vidard) Initial version !!---------------------------------------------------------------------- !! * Arguments REAL(wp), INTENT( in ) :: & ptem, & ! potential temperature psal ! salinity REAL(wp), INTENT( in ) :: & ptem_tl, & ! potential temperature psal_tl ! salinity REAL(wp), INTENT( out ) :: & prhop_tl ! potential density (surface referenced) !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & ! temporary scalars & zt, zs, zsr, zr1, zr2, zr3, zr4, zrhop, & & zttl, zstl, zsrtl, zr1tl, zr2tl, zr3tl, zr4tl, zrhoptl REAL(wp) :: zws, zwstl, zeps !!---------------------------------------------------------------------- #ifdef key_sp zeps = 1.e-7 #else zeps = 1.e-14 #endif SELECT CASE ( nn_eos ) CASE ( 0 ) ! Jackett and McDougall (1994) formulation zws = SQRT( ABS( psal ) ) zt = ptem zs = psal ! square root salinity zsr= zws ! compute volumic mass pure water at atm pressure zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 ! seawater volumic mass atm pressure zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt & -4.0899e-3 ) *zt+0.824493 zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 zr4= 4.8314e-4 ! ================== ! tangent linear part ! ================== zwstl = 1/MAX( 2*zws , zeps ) * psal_tl zttl = ptem_tl zstl = psal_tl ! square root salinity zsrtl= zwstl ! compute volumic mass pure water at atm pressure zr1tl= ( ( ( ( 5.*6.536332e-9 * zt & & -4.*1.120083e-6 ) * zt & & +3.*1.001685e-4 ) * zt & & -2.*9.095290e-3 ) * zt & & + 6.793952e-2 ) * zttl ! seawater volumic mass atm pressure zr2tl= ( ( ( 4.*5.3875e-9 * zt & & -3.*8.2467e-7 ) * zt & & +2.*7.6438e-5 ) * zt & & - 4.0899e-3 ) * zttl zr3tl= ( -2.*1.6546e-6 * zt & & + 1.0227e-4 ) * zttl ! potential volumic mass (reference to the surface) zrhoptl= zr1tl & & + zs * zr2tl & & + zsr * zs * zr3tl & & + zr3 * zs * zsrtl & & + ( 2. * zr4 * zs + zr2 & & + zr3 * zsr ) * zstl ! save potential volumic mass prhop_tl = zrhoptl CASE ( 1 ) ! Linear formulation function of temperature only zttl = ptem_tl ! ... potential volumic mass prhop_tl = ( rau0 * ( - rn_alpha * zttl ) ) CASE ( 2 ) ! Linear formulation function of temperature and salinity zttl = ptem_tl zstl = psal_tl ! ... potential volumic mass prhop_tl = rau0 * ( rn_beta * zstl - rn_alpha * zttl ) CASE DEFAULT WRITE(ctmp1,*) ' bad flag value for nn_eos = ', nn_eos CALL ctl_stop( ctmp1 ) END SELECT END SUBROUTINE eos_pot_1pt_tan SUBROUTINE eos_insitu_2d_tan( ptem, psal, pdep, ptem_tl, psal_tl, prd_tl ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_insitu_2d_tan : TL OF ROUTINE eos_insitu_2d *** !! !! ** Purpose of direct routine : Compute the in situ density !! (ratio rho/rau0) from potential temperature and salinity !! using an equation of state defined through the namelist !! parameter nn_eos. * 2D field case !! !! ** Method of direct routine : 3 cases: !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. !! the in situ density is computed directly as a function of !! potential temperature relative to the surface (the opa t !! variable), salt and pressure (assuming no pressure variation !! along geopotential surfaces, i.e. the pressure p in decibars !! is approximated by the depth in meters. !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 !! with pressure p decibars !! potential temperature t deg celsius !! salinity s psu !! reference volumic mass rau0 kg/m**3 !! in situ volumic mass rho kg/m**3 !! in situ density anomalie prd no units !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, !! t = 40 deg celcius, s=40 psu !! nn_eos = 1 : linear equation of state function of temperature only !! prd(t) = 0.0285 - ralpha * t !! nn_eos = 2 : linear equation of state function of temperature and !! salinity !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. !! Note that no boundary condition problem occurs in this routine !! as (ptem,psal) are defined over the whole domain. !! !! ** Comments on Adjoint Routine : !! Care has been taken to avoid division by zero when computing !! the inverse of the square root of salinity at masked salinity !! points. !! !! ** Action : !! !! References : !! !! History : !! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eostan.F !! 9.0 ! 07-07 (K. Mogensen) Initial version based on eostan.F !! ! 08-07 (A. Vidard) bug fix in computation of prd_tl if neos=1 !!----------------------------------------------------------------------- !! * Modules used !! * Arguments REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: & & ptem, & ! potential temperature & psal, & ! salinity & pdep, & ! depth & ptem_tl, & ! TL of potential temperature & psal_tl ! TL of salinity REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) :: & & prd_tl ! TL of potential density (surface referenced) !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & ! temporary scalars zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, & zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, & zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, & zr4tl, zrhoptl, zetl, zbwtl, & zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, & zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, & zmask REAL(wp), DIMENSION(jpi,jpj) :: zws !!---------------------------------------------------------------------- SELECT CASE ( nn_eos ) CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==! #ifdef key_sp zeps = 1.e-7 #else zeps = 1.e-14 #endif !CDIR NOVERRCHK DO jj = 1, jpjm1 !CDIR NOVERRCHK DO ji = 1, fs_jpim1 ! vector opt. zws(ji,jj) = SQRT( ABS( psal(ji,jj) ) ) END DO END DO DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask zt = ptem (ji,jj) ! interpolated T zs = psal (ji,jj) ! interpolated S zsr= zws(ji,jj) ! square root of interpolated S zh = pdep(ji,jj) ! depth at the partial step level ! compute volumic mass pure water at atm pressure zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt & & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 ! seawater volumic mass atm pressure zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt & & -4.0899e-3 ) *zt+0.824493 zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 zr4= 4.8314e-4 ! potential volumic mass (reference to the surface) zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 ! add the compression terms ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 zbw= ( 1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 zb = zbw + ze * zs zd = -2.042967e-2 zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 za = ( zd*zsr + zc ) *zs + zaw zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545 za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 zk0= ( zb1*zsr + za1 )*zs + zkw ! Tangent linear part zttl = ptem_tl(ji,jj) zstl = psal_tl(ji,jj) zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) & & * tmask(ji,jj,1) * zstl zr1tl= ( ( ( ( 5.*6.536332e-9 * zt & & -4.*1.120083e-6 ) * zt & & +3.*1.001685e-4 ) * zt & & -2.*9.095290e-3 ) * zt & & + 6.793952e-2 ) * zttl zr2tl= ( ( ( 4.*5.3875e-9 * zt & & -3.*8.2467e-7 ) * zt & & +2.*7.6438e-5 ) * zt & & - 4.0899e-3 ) * zttl zr3tl= ( -2.*1.6546e-6 * zt & & + 1.0227e-4 ) * zttl zrhoptl= zr1tl & & + zs * zr2tl & & + zsr * zs * zr3tl & & + zr3 * zs * zsrtl & & + ( 2. * zr4 * zs + zr2 & & + zr3 * zsr ) * zstl zetl = ( -2.*3.508914e-8 * zt & & - 1.248266e-8 ) * zttl zbwtl= ( 2.*1.296821e-6 * zt & & - 5.782165e-9 ) * zttl zbtl = zbwtl & & + zs * zetl & & + ze * zstl zctl = ( -2.*7.267926e-5 * zt & & + 2.598241e-3 ) * zttl zawtl= ( ( 3.*5.939910e-6 * zt & & +2.*2.512549e-3 ) * zt & & - 0.1028859 ) * zttl zatl = zawtl & & + zd * zs * zsrtl & & + zs * zctl & & + ( zd * zsr + zc ) * zstl zb1tl= ( -2.*0.1909078 * zt & & + 7.390729 ) * zttl za1tl= ( ( 3.*2.326469e-3 * zt & & +2.*1.553190 ) * zt & & - 65.00517 ) * zttl zkwtl= ( ( ( -4.*1.361629e-4 * zt & & -3.*1.852732e-2 ) * zt & & -2.*30.41638 ) * zt & & + 2098.925 ) * zttl zk0tl= zkwtl & & + zb1 * zs * zsrtl & & + zs * zsr * zb1tl & & + zs * za1tl & & + ( zb1 * zsr + za1 ) * zstl ! Masked in situ density anomaly zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) ) zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 ) prd_tl(ji,jj) = tmask(ji,jj,1) * zrdc2 * & & ( zrhoptl & & - zrdc2 * zh * zrdc1**2 * zrhop & & * ( zk0tl & & - zh * ( zatl & & - zh * zbtl ) ) )& & / rau0 END DO END DO ! CASE ( 1 ) !== Linear formulation = F( temperature ) ==! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. prd_tl(ji,jj) = ( - rn_alpha * ptem_tl(ji,jj) ) * tmask(ji,jj,1) END DO END DO ! CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. prd_tl (ji,jj) = ( rn_beta * psal_tl(ji,jj) - rn_alpha * ptem_tl(ji,jj) ) * tmask(ji,jj,1) END DO END DO ! END SELECT END SUBROUTINE eos_insitu_2d_tan SUBROUTINE eos_insitu_adj(ptem, psal, ptem_ad, psal_ad, prd_ad) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_insitu_tan : Adjoint OF ROUTINE eos_insitu *** !! !! ** Purpose of direct routine : Compute the in situ density !! (ratio rho/rau0) from potential temperature and salinity !! using an equation of state defined through the namelist !! parameter nneos. !! !! ** Method of direct routine : 3 cases: !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. !! the in situ density is computed directly as a function of !! potential temperature relative to the surface (the opa t !! variable), salt and pressure (assuming no pressure variation !! along geopotential surfaces, i.e. the pressure p in decibars !! is approximated by the depth in meters. !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 !! with pressure p decibars !! potential temperature t deg celsius !! salinity s psu !! reference volumic mass rau0 kg/m**3 !! in situ volumic mass rho kg/m**3 !! in situ density anomalie prd no units !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, !! t = 40 deg celcius, s=40 psu !! nn_eos = 1 : linear equation of state function of temperature only !! prd(t) = 0.0285 - rn_alpha * t !! nn_eos = 2 : linear equation of state function of temperature and !! salinity !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. !! Note that no boundary condition problem occurs in this routine !! as (ptem,psal) are defined over the whole domain. !! !! ** Comments on Adjoint Routine : !! Care has been taken to avoid division by zero when computing !! the inverse of the square root of salinity at masked salinity !! points. !! !! ** Action : !! !! References : !! !! History : !! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eostan.F !! 9.0 ! 08-08 (A. Vidard) 9.0 version !!----------------------------------------------------------------------- !! * Modules used !! * Arguments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & ptem, & ! potential temperature psal ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & ptem_ad, & ! potential temperature psal_ad ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & prd_ad ! potential density (surface referenced) !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & ! temporary scalars zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, & zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, & ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, & zr4ad, zrhopad, zead, zbwad, & zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, & zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, & zmask, zrau0r REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws !!---------------------------------------------------------------------- ! initialization of adjoint variables ztad = 0.0_wp zsad = 0.0_wp zhad = 0.0_wp zsrad = 0.0_wp zr1ad = 0.0_wp zr2ad = 0.0_wp zr3ad = 0.0_wp zr4ad = 0.0_wp zrhopad = 0.0_wp zead = 0.0_wp zbwad = 0.0_wp zbad = 0.0_wp zdad = 0.0_wp zcad = 0.0_wp zawad = 0.0_wp zaad = 0.0_wp zb1ad = 0.0_wp za1ad = 0.0_wp zkwad = 0.0_wp zk0ad = 0.0_wp SELECT CASE ( nn_eos ) CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==! zrau0r = 1.e0 / rau0 #ifdef key_sp zeps = 1.e-7 #else zeps = 1.e-14 #endif !CDIR NOVERRCHK zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) ) DO jk = 1, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zt = ptem(ji,jj,jk) zs = psal(ji,jj,jk) zh = fsdept(ji,jj,jk) ! depth zsr= zws(ji,jj,jk) ! square root salinity ! compute volumic mass pure water at atm pressure zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 ! seawater volumic mass atm pressure zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt & -4.0899e-3 ) *zt+0.824493 zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 zr4= 4.8314e-4 ! potential volumic mass (reference to the surface) zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 ! add the compression terms ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 zbw= ( 1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 zb = zbw + ze * zs zd = -2.042967e-2 zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 za = ( zd*zsr + zc ) *zs + zaw zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545 za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 zk0= ( zb1*zsr + za1 )*zs + zkw zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) ) zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 ) ! ============ ! Adjoint part ! ============ ! Masked in situ density anomaly zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) & & * zrdc2 * zrau0r zk0ad = zk0ad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & * zrau0r zaad = zaad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & * zh * zrau0r zbad = zbad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & * zh * zh * zrau0r prd_ad(ji,jj,jk) = 0.0_wp zkwad = zkwad + zk0ad zsrad = zsrad + zk0ad * zb1 * zs zb1ad = zb1ad + zk0ad * zs * zsr za1ad = za1ad + zk0ad * zs zsad = zsad + zk0ad * ( zb1 * zsr + za1 ) zk0ad = 0.0_wp ztad = ztad + zkwad * ( ( (-4.*1.361629e-4 * zt & & -3.*1.852732e-2 ) * zt & & -2.*30.41638 ) * zt & & + 2098.925 ) zkwad = 0.0_wp ztad = ztad + za1ad * ( ( 3.*2.326469e-3 * zt & & +2.*1.553190 ) * zt & & - 65.00517 ) za1ad = 0.0_wp ztad = ztad + zb1ad * (-2.*0.1909078 * zt & & + 7.390729 ) zb1ad = 0.0_wp zawad = zawad + zaad zsrad = zsrad + zaad * zd * zs zcad = zcad + zaad * zs zsad = zsad + zaad * ( zd * zsr + zc ) zaad = 0.0_wp ztad = ztad + zawad * ( ( 3.*5.939910e-6 * zt & & +2.*2.512549e-3 ) * zt & & - 0.1028859 ) zawad = 0.0_wp ztad = ztad + zcad * (-2.*7.267926e-5 * zt & & + 2.598241e-3 ) zcad = 0.0_wp zbwad = zbwad + zbad zead = zead + zbad * zs zsad = zsad + zbad * ze zbad = 0.0_wp ztad = ztad + zbwad * ( 2.*1.296821e-6 * zt & & - 5.782165e-9 ) zbwad = 0.0_wp ztad = ztad + zead * (-2.*3.508914e-8 * zt & & - 1.248266e-8 ) zead = 0.0_wp zr1ad = zr1ad + zrhopad zr2ad = zr2ad + zrhopad * zs zr3ad = zr3ad + zrhopad * zsr * zs zsrad = zsrad + zrhopad * zr3 * zs zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 & & + zr3 * zsr ) zrhopad = 0.0_wp ztad = ztad + zr3ad * (-2.*1.6546e-6 * zt & & + 1.0227e-4 ) zr3ad = 0.0_wp ztad = ztad + zr2ad * ( ( ( 4.*5.3875e-9 * zt & & -3.*8.2467e-7 ) * zt & & +2.*7.6438e-5 ) * zt & & - 4.0899e-3 ) zr2ad = 0.0_wp ztad = ztad + zr1ad * ( ( ( ( 5.*6.536332e-9 * zt & & -4.*1.120083e-6 ) * zt & & +3.*1.001685e-4 ) * zt & & -2.*9.095290e-3 ) * zt & & + 6.793952e-2 ) zr1ad = 0.0_wp zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) & & * tmask(ji,jj,jk) zsrad = 0.0_wp psal_ad(ji,jj,jk) = psal_ad(ji,jj,jk) + zsad ptem_ad(ji,jj,jk) = ptem_ad(ji,jj,jk) + ztad ztad = 0.0_wp zsad = 0.0_wp END DO END DO END DO ! CASE ( 1 ) !== Linear formulation function of temperature only ==! DO jk = 1, jpkm1 ptem_ad(:,:,jk) = ptem_ad(:,:,jk) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk) prd_ad(:,:,jk) = 0.0_wp END DO ! CASE ( 2 ) !== Linear formulation function of temperature and salinity ==! DO jk = 1, jpkm1 ptem_ad(:,:,jk) = ptem_ad(:,:,jk) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk) psal_ad(:,:,jk) = psal_ad(:,:,jk) + rn_beta * prd_ad( :,:,jk) * tmask(:,:,jk) prd_ad( :,:,jk) = 0.0_wp END DO ! END SELECT END SUBROUTINE eos_insitu_adj SUBROUTINE eos_insitu_pot_adj ( ptem, psal, ptem_ad, psal_ad, prd_ad, prhop_ad ) !!---------------------------------------------------------------------- !! *** ROUTINE eos_insitu_pot_adj *** !! !! ** Purpose or the direct routine: !! Compute the in situ density (ratio rho/rau0) and the !! potential volumic mass (Kg/m3) from potential temperature and !! salinity fields using an equation of state defined through the !! namelist parameter nn_eos. !! !! ** Method : !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. !! the in situ density is computed directly as a function of !! potential temperature relative to the surface (the opa t !! variable), salt and pressure (assuming no pressure variation !! along geopotential surfaces, i.e. the pressure p in decibars !! is approximated by the depth in meters. !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 !! rhop(t,s) = rho(t,s,0) !! with pressure p decibars !! potential temperature t deg celsius !! salinity s psu !! reference volumic mass rau0 kg/m**3 !! in situ volumic mass rho kg/m**3 !! in situ density anomalie prd no units !! !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, !! t = 40 deg celcius, s=40 psu !! !! neos = 1 : linear equation of state function of temperature only !! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - ralpha * t !! rhop(t,s) = rho(t,s) !! !! nn_eos = 2 : linear equation of state function of temperature and !! salinity !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 !! = rn_beta * s - rn_alpha * tn - 1. !! rhop(t,s) = rho(t,s) !! Note that no boundary condition problem occurs in this routine !! as (tn,sn) or (ta,sa) are defined over the whole domain. !! !! ** Action : - prd , the in situ density (no units) !! - prhop, the potential volumic mass (Kg/m3) !! !! References : !! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994 !! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978 !! !! History of the adjoint routine: !! 9.0 ! 08-06 (A. Vidard) Initial version !!---------------------------------------------------------------------- !! * Arguments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & ptem, & ! potential temperature psal ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & ptem_ad, & ! potential temperature psal_ad ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & prd_ad, & ! potential density (surface referenced) prhop_ad !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & ! temporary scalars zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, & zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, & ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, & zr4ad, zrhopad, zead, zbwad, & zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, & zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, & zmask, zrau0r REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws !!---------------------------------------------------------------------- ! initialization of adjoint variables ztad = 0.0_wp zsad = 0.0_wp zhad = 0.0_wp zsrad = 0.0_wp zr1ad = 0.0_wp zr2ad = 0.0_wp zr3ad = 0.0_wp zr4ad = 0.0_wp zrhopad = 0.0_wp zead = 0.0_wp zbwad = 0.0_wp zbad = 0.0_wp zdad = 0.0_wp zcad = 0.0_wp zawad = 0.0_wp zaad = 0.0_wp zb1ad = 0.0_wp za1ad = 0.0_wp zkwad = 0.0_wp zk0ad = 0.0_wp SELECT CASE ( nn_eos ) CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==! zrau0r = 1.e0 / rau0 #ifdef key_sp zeps = 1.e-7 #else zeps = 1.e-14 #endif !CDIR NOVERRCHK zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) ) ! DO jk = jpkm1, 1, -1 DO jj = jpj, 1, -1 DO ji = jpi, 1, -1 ! direct recomputing zt = ptem(ji,jj,jk) zs = psal(ji,jj,jk) zh = fsdept(ji,jj,jk) ! depth zsr = zws(ji,jj,jk) ! square root salinity ! compute volumic mass pure water at atm pressure zr1 = ( ( ( ( 6.536332e-9 * zt - 1.120083e-6 ) * zt & & + 1.001685e-4 ) * zt - 9.095290e-3 ) * zt & & + 6.793952e-2 ) * zt + 999.842594 ! seawater volumic mass atm pressure zr2 = ( ( ( 5.3875e-9 * zt - 8.2467e-7 ) * zt & & + 7.6438e-5 ) * zt - 4.0899e-3 ) * zt + 0.824493 zr3 = ( -1.6546e-6 * zt + 1.0227e-4 ) * zt - 5.72466e-3 zr4 = 4.8314e-4 ! potential volumic mass (reference to the surface) zrhop = ( zr4 * zs + zr3*zsr + zr2 ) * zs + zr1 ! add the compression terms ze = ( -3.508914e-8 * zt - 1.248266e-8 ) * zt - 2.595994e-6 zbw = ( 1.296821e-6 * zt - 5.782165e-9 ) * zt + 1.045941e-4 zb = zbw + ze * zs zd = -2.042967e-2 zc = (-7.267926e-5 * zt + 2.598241e-3 ) * zt + 0.1571896 zaw= ( ( 5.939910e-6 * zt + 2.512549e-3 ) * zt - 0.1028859 & & ) * zt - 4.721788 za = ( zd * zsr + zc ) * zs + zaw zb1= (-0.1909078 * zt + 7.390729 ) * zt - 55.87545 za1= ( ( 2.326469e-3 * zt + 1.553190 ) * zt - 65.00517 & & ) * zt + 1044.077 zkw= ( ( (-1.361629e-4 * zt - 1.852732e-2 ) * zt - 30.41638 & & ) * zt + 2098.925 ) * zt + 190925.6 zk0= ( zb1 * zsr + za1 ) * zs + zkw zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) ) zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 ) ! ============ ! Adjoint part ! ============ ! Masked in situ density anomaly zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) & & * zrdc2 * zrau0r zk0ad = zk0ad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & * zrau0r zaad = zaad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & * zh * zrau0r zbad = zbad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & * zh * zh * zrau0r prd_ad(ji,jj,jk) = 0.0_wp zkwad = zkwad + zk0ad zsrad = zsrad + zk0ad * zb1 * zs zb1ad = zb1ad + zk0ad * zs * zsr za1ad = za1ad + zk0ad * zs zsad = zsad + zk0ad * ( zb1 * zsr + za1 ) zk0ad = 0.0_wp ztad = ztad + zkwad * ( ( (-4.*1.361629e-4 * zt & & -3.*1.852732e-2 ) * zt & & -2.*30.41638 ) * zt & & + 2098.925 ) zkwad = 0.0_wp ztad = ztad + za1ad * ( ( 3.*2.326469e-3 * zt & & +2.*1.553190 ) * zt & & - 65.00517 ) za1ad = 0.0_wp ztad = ztad + zb1ad * (-2.*0.1909078 * zt & & + 7.390729 ) zb1ad = 0.0_wp zawad = zawad + zaad zsrad = zsrad + zaad * zd * zs zcad = zcad + zaad * zs zsad = zsad + zaad * ( zd * zsr + zc ) zaad = 0.0_wp ztad = ztad + zawad * ( ( 3.*5.939910e-6 * zt & & +2.*2.512549e-3 ) * zt & & - 0.1028859 ) zawad = 0.0_wp ztad = ztad + zcad * (-2.*7.267926e-5 * zt & & + 2.598241e-3 ) zcad = 0.0_wp zsad = zsad + zbad * ze zead = zead + zbad * zs zbwad = zbwad + zbad zbad = 0.0_wp ztad = ztad + zbwad * ( 2.*1.296821e-6 * zt & & - 5.782165e-9 ) zbwad = 0.0_wp ztad = ztad + zead * (-2.*3.508914e-8 * zt & & - 1.248266e-8 ) zead = 0.0_wp zrhopad = zrhopad + prhop_ad(ji,jj,jk) * tmask(ji,jj,jk) prhop_ad(ji,jj,jk) = 0.0_wp zr1ad = zr1ad + zrhopad zr2ad = zr2ad + zrhopad * zs zr3ad = zr3ad + zrhopad * zsr * zs zsrad = zsrad + zrhopad * zr3 * zs zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 & & + zr3 * zsr ) zrhopad = 0.0_wp ztad = ztad + zr3ad * (-2.*1.6546e-6 * zt & & + 1.0227e-4 ) zr3ad = 0.0_wp ztad = ztad + zr2ad * ( ( ( 4.*5.3875e-9 * zt & & -3.*8.2467e-7 ) * zt & & +2.*7.6438e-5 ) * zt & & - 4.0899e-3 ) zr2ad = 0.0_wp ztad = ztad + zr1ad * ( ( ( ( 5.*6.536332e-9 * zt & & -4.*1.120083e-6 ) * zt & & +3.*1.001685e-4 ) * zt & & -2.*9.095290e-3 ) * zt & & + 6.793952e-2 ) zr1ad = 0.0_wp zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) & & * tmask(ji,jj,jk) zsrad = 0.0_wp psal_ad(ji,jj,jk) = psal_ad(ji,jj,jk) + zsad ptem_ad(ji,jj,jk) = ptem_ad(ji,jj,jk) + ztad ztad = 0.0_wp zsad = 0.0_wp END DO END DO END DO ! CASE ( 1 ) !== Linear formulation = F( temperature ) ==! DO jk = jpkm1, 1, -1 prd_ad(:,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk) prhop_ad(:,:,jk) = 0.0_wp ptem_ad(:,:,jk) = ptem_ad(:,:,jk) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk) prd_ad(:,:,jk) = 0.0_wp END DO ! CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! DO jk = 1, jpkm1 prd_ad( :,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk) prhop_ad(:,:,jk) = 0.0_wp ptem_ad( :,:,jk) = ptem_ad(:,:,jk) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk) psal_ad( :,:,jk) = psal_ad(:,:,jk) + rn_beta * prd_ad(:,:,jk) * tmask(:,:,jk) prd_ad( :,:,jk) = 0.0_wp END DO ! END SELECT END SUBROUTINE eos_insitu_pot_adj SUBROUTINE eos_pot_1pt_adj( ptem, psal, ptem_ad, psal_ad, prhop_ad ) REAL(wp), INTENT( in ) :: & ptem, & ! potential temperature psal ! salinity REAL(wp), INTENT( inout ) :: & ptem_ad, & ! potential temperature psal_ad ! salinity REAL(wp), INTENT( inout ) :: & prhop_ad ! potential density (surface referenced) !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & ! temporary scalars & zt, zs, zsr, zr1, zr2, zr3, zr4, zrhop, & & ztad, zsad, zsrad, zr1ad, zr2ad, zr3ad, zr4ad, zrhopad REAL(wp) :: zws, zwsad, zeps !!---------------------------------------------------------------------- #ifdef key_sp zeps = 1.e-7 #else zeps = 1.e-14 #endif zwsad = 0.0_wp ztad = 0.0_wp zsad = 0.0_wp zsrad = 0.0_wp zr1ad = 0.0_wp zr2ad = 0.0_wp zr3ad = 0.0_wp zr4ad = 0.0_wp zrhopad = 0.0_wp SELECT CASE ( nn_eos ) CASE ( 0 ) ! Jackett and McDougall (1994) formulation zws = SQRT( ABS( psal ) ) zt = ptem zs = psal ! square root salinity zsr= zws ! compute volumic mass pure water at atm pressure zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 ! seawater volumic mass atm pressure zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt & -4.0899e-3 ) *zt+0.824493 zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 zr4= 4.8314e-4 ! ============ ! Adjoint part ! ============ ! save potential volumic mass zrhopad = zrhopad + prhop_ad prhop_ad = 0.0_wp ! potential volumic mass (reference to the surface) zr1ad = zr1ad + zrhopad zr2ad = zr2ad + zrhopad * zs zr3ad = zr3ad + zrhopad * zsr * zs zsrad = zsrad + zrhopad * zr3 * zs zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 & & + zr3 * zsr ) zrhopad= 0.0_wp ! seawater volumic mass atm pressure ztad = ztad + zr3ad * ( -2.*1.6546e-6 * zt & & + 1.0227e-4 ) zr3ad = 0.0_wp ztad = ztad + zr2ad * ( ( ( 4.*5.3875e-9 * zt & & -3.*8.2467e-7 ) * zt & & +2.*7.6438e-5 ) * zt & & - 4.0899e-3 ) zr2ad = 0.0_wp ! compute volumic mass pure water at atm pressure ztad = ztad + zr1ad * ( ( ( ( 5.*6.536332e-9 * zt & & -4.*1.120083e-6 ) * zt & & +3.*1.001685e-4 ) * zt & & -2.*9.095290e-3 ) * zt & & + 6.793952e-2 ) zr1ad = 0.0_wp ! square root salinity zwsad = zwsad + zsrad zsrad = 0.0_wp ptem_ad = ptem_ad + ztad psal_ad = psal_ad + zsad ztad = 0.0_wp zsad = 0.0_wp psal_ad = psal_ad + 1/MAX( 2*zws , zeps ) * zwsad zwsad = 0.0_wp CASE ( 1 ) ! Linear formulation function of temperature only ! ... potential volumic mass ptem_ad = ptem_ad - prhop_ad * rau0 * rn_alpha prhop_ad = 0.0_wp CASE ( 2 ) ! Linear formulation function of temperature and salinity ! ... potential volumic mass ptem_ad = ptem_ad - prhop_ad * rau0 * rn_alpha psal_ad = psal_ad + prhop_ad * rau0 * rn_beta prhop_ad = 0.0_wp CASE DEFAULT WRITE(ctmp1,*) ' bad flag value for nn_eos = ', nn_eos CALL ctl_stop( ctmp1 ) END SELECT END SUBROUTINE eos_pot_1pt_adj SUBROUTINE eos_insitu_2d_adj( ptem, psal, pdep, ptem_ad, psal_ad, prd_ad ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_insitu_2d_adj : adj OF ROUTINE eos_insitu_2d *** !! !! ** Purpose of direct routine : Compute the in situ density !! (ratio rho/rau0) from potential temperature and salinity !! using an equation of state defined through the namelist !! parameter nn_eos. * 2D field case !! !! ** Method of direct routine : 3 cases: !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. !! the in situ density is computed directly as a function of !! potential temperature relative to the surface (the opa t !! variable), salt and pressure (assuming no pressure variation !! along geopotential surfaces, i.e. the pressure p in decibars !! is approximated by the depth in meters. !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 !! with pressure p decibars !! potential temperature t deg celsius !! salinity s psu !! reference volumic mass rau0 kg/m**3 !! in situ volumic mass rho kg/m**3 !! in situ density anomalie prd no units !! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar, !! t = 40 deg celcius, s=40 psu !! nn_eos = 1 : linear equation of state function of temperature only !! prd(t) = 0.0285 - rn_alpha * t !! nn_eos = 2 : linear equation of state function of temperature and !! salinity !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. !! Note that no boundary condition problem occurs in this routine !! as (ptem,psal) are defined over the whole domain. !! !! ** Comments on Adjoint Routine : !! Care has been taken to avoid division by zero when computing !! the inverse of the square root of salinity at masked salinity !! points. !! !! ** Action : !! !! References : !! !! History : !! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eosadj.F !! 9.0 ! 08-07 (A. Vidard) first version based on eosadj !!----------------------------------------------------------------------- !! * Modules used !! * Arguments REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: & & ptem, & ! potential temperature & psal, & ! salinity & pdep ! depth REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: & & ptem_ad, & ! TL of potential temperature & psal_ad ! TL of salinity REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: & & prd_ad ! TL of potential density (surface referenced) INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: & ! temporary scalars zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, & zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, & ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, & zr4ad, zrhopad, zead, zbwad, & zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, & zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, & zmask REAL(wp), DIMENSION(jpi,jpj) :: zws !!---------------------------------------------------------------------- ! initialization of adjoint variables ztad = 0.0_wp zsad = 0.0_wp zhad = 0.0_wp zsrad = 0.0_wp zr1ad = 0.0_wp zr2ad = 0.0_wp zr3ad = 0.0_wp zr4ad = 0.0_wp zrhopad = 0.0_wp zead = 0.0_wp zbwad = 0.0_wp zbad = 0.0_wp zdad = 0.0_wp zcad = 0.0_wp zawad = 0.0_wp zaad = 0.0_wp zb1ad = 0.0_wp za1ad = 0.0_wp zkwad = 0.0_wp zk0ad = 0.0_wp SELECT CASE ( nn_eos ) CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==! #ifdef key_sp zeps = 1.e-7 #else zeps = 1.e-14 #endif !CDIR NOVERRCHK DO jj = 1, jpjm1 !CDIR NOVERRCHK DO ji = 1, fs_jpim1 ! vector opt. zws(ji,jj) = SQRT( ABS( psal(ji,jj) ) ) END DO END DO ! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask zt = ptem (ji,jj) ! interpolated T zs = psal (ji,jj) ! interpolated S zsr= zws(ji,jj) ! square root of interpolated S zh = pdep(ji,jj) ! depth at the partial step level ! compute volumic mass pure water at atm pressure zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt & & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 ! seawater volumic mass atm pressure zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt & & -4.0899e-3 ) *zt+0.824493 zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 zr4= 4.8314e-4 ! potential volumic mass (reference to the surface) zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 ! add the compression terms ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 zbw= ( 1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 zb = zbw + ze * zs zd = -2.042967e-2 zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 za = ( zd*zsr + zc ) *zs + zaw zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545 za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 zk0= ( zb1*zsr + za1 )*zs + zkw zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) ) zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 ) ! ============ ! Adjoint part ! ============ ! Masked in situ density anomaly zrhopad = zrhopad + prd_ad(ji,jj) * tmask(ji,jj,1) & & * zrdc2 / rau0 zk0ad = zk0ad - prd_ad(ji,jj) * tmask(ji,jj,1) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & / rau0 zaad = zaad + prd_ad(ji,jj) * tmask(ji,jj,1) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & * zh / rau0 zbad = zbad - prd_ad(ji,jj) * tmask(ji,jj,1) & & * zrdc2 * zrdc2 * zh & & * zrdc1**2 * zrhop & & * zh * zh / rau0 prd_ad(ji,jj) = 0.0_wp zkwad = zkwad + zk0ad zsrad = zsrad + zk0ad * zb1 * zs zb1ad = zb1ad + zk0ad * zs * zsr za1ad = za1ad + zk0ad * zs zsad = zsad + zk0ad * ( zb1 * zsr + za1 ) zk0ad = 0.0_wp ztad = ztad + zkwad * ( ( (-4.*1.361629e-4 * zt & & -3.*1.852732e-2 ) * zt & & -2.*30.41638 ) * zt & & + 2098.925 ) zkwad = 0.0_wp ztad = ztad + za1ad * ( ( 3.*2.326469e-3 * zt & & +2.*1.553190 ) * zt & & - 65.00517 ) za1ad = 0.0_wp ztad = ztad + zb1ad * (-2.*0.1909078 * zt & & + 7.390729 ) zb1ad = 0.0_wp zawad = zawad + zaad zsrad = zsrad + zaad * zd * zs zcad = zcad + zaad * zs zsad = zsad + zaad * ( zd * zsr + zc ) zaad = 0.0_wp ztad = ztad + zawad * ( ( 3.*5.939910e-6 * zt & & +2.*2.512549e-3 ) * zt & & - 0.1028859 ) zawad = 0.0_wp ztad = ztad + zcad * (-2.*7.267926e-5 * zt & & + 2.598241e-3 ) zcad = 0.0_wp zbwad = zbwad + zbad zead = zead + zbad * zs zsad = zsad + zbad * ze zbad = 0.0_wp ztad = ztad + zbwad * ( 2.*1.296821e-6 * zt & & - 5.782165e-9 ) zbwad = 0.0_wp ztad = ztad + zead * (-2.*3.508914e-8 * zt & & - 1.248266e-8 ) zead = 0.0_wp zr1ad = zr1ad + zrhopad zr2ad = zr2ad + zrhopad * zs zr3ad = zr3ad + zrhopad * zsr * zs zsrad = zsrad + zrhopad * zr3 * zs zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 & & + zr3 * zsr ) zrhopad = 0.0_wp ztad = ztad + zr3ad * (-2.*1.6546e-6 * zt & & + 1.0227e-4 ) zr3ad = 0.0_wp ztad = ztad + zr2ad * ( ( ( 4.*5.3875e-9 * zt & & -3.*8.2467e-7 ) * zt & & +2.*7.6438e-5 ) * zt & & - 4.0899e-3 ) zr2ad = 0.0_wp ztad = ztad + zr1ad * ( ( ( ( 5.*6.536332e-9 * zt & & -4.*1.120083e-6 ) * zt & & +3.*1.001685e-4 ) * zt & & -2.*9.095290e-3 ) * zt & & + 6.793952e-2 ) zr1ad = 0.0_wp zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) & & * tmask(ji,jj, 1) zsrad = 0.0_wp psal_ad(ji,jj) = psal_ad(ji,jj) + zsad ptem_ad(ji,jj) = ptem_ad(ji,jj) + ztad ztad = 0.0_wp zsad = 0.0_wp END DO END DO ! CASE ( 1 ) !== Linear formulation = F( temperature ) ==! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ptem_ad(ji,jj) = ptem_ad(ji,jj) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1) prd_ad(ji,jj) = 0.0_wp END DO END DO ! CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ptem_ad(ji,jj) = ptem_ad(ji,jj) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1) psal_ad(ji,jj) = psal_ad(ji,jj) + prd_ad(ji,jj) * rn_beta * tmask(ji,jj,1) prd_ad (ji,jj) = 0.0_wp END DO END DO ! END SELECT END SUBROUTINE eos_insitu_2d_adj SUBROUTINE eos_bn2_tan ( ptem, psal, ptem_tl, psal_tl, pn2_tl ) !!---------------------------------------------------------------------- !! *** ROUTINE eos_bn2_tan *** !! !! ** Purpose of the direct routine: Compute the local !! Brunt-Vaisala frequency at the time-step of the input arguments !! !! ** Method of the direct routine: !! * nn_eos = 0 : UNESCO sea water properties !! The brunt-vaisala frequency is computed using the polynomial !! polynomial expression of McDougall (1987): !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is !! computed and used in zdfddm module : !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) !! * nn_eos = 1 : linear equation of state (temperature only) !! N^2 = grav * rn_alpha * dk[ t ]/e3w !! * nn_eos = 2 : linear equation of state (temperature & salinity) !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w !! The use of potential density to compute N^2 introduces e r r o r !! in the sign of N^2 at great depths. We recommand the use of !! nn_eos = 0, except for academical studies. !! Macro-tasked on horizontal slab (jk-loop) !! N.B. N^2 is set to zero at the first level (JK=1) in inidtr !! and is never used at this level. !! !! ** Action : - pn2 : the brunt-vaisala frequency !! !! References : !! McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987. !! !! History: !! ! 08-07 (A. Vidard) First version !!---------------------------------------------------------------------- !! * Arguments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & ptem, & ! potential temperature psal ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & ptem_tl, & ! potential temperature psal_tl ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: & pn2_tl ! Brunt-Vaisala frequency !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zgde3w, zt, zs, zh, & ! temporary scalars zalbet, zbeta ! " " REAL(wp) :: & zttl, zstl, & ! temporary scalars zalbettl, zbetatl ! " " #if defined key_zdfddm REAL(wp) :: zds, zdstl ! temporary scalars #endif ! pn2_tl : interior points only (2=< jk =< jpkm1 ) ! -------------------------- SELECT CASE ( nn_eos ) CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==! DO jk = 2, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zgde3w = grav / fse3w(ji,jj,jk) zt = 0.5 * ( ptem(ji,jj,jk) + ptem(ji,jj,jk-1) ) ! potential temperature at w-point zs = 0.5 * ( psal(ji,jj,jk) + psal(ji,jj,jk-1) ) - 35.0 ! salinity anomaly (s-35) at w-point zh = fsdepw(ji,jj,jk) ! depth in meters at w-point zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta & - 0.203814e-03 ) * zt & & + 0.170907e-01 ) * zt & & + 0.665157e-01 & & + ( - 0.678662e-05 * zs & & - 0.846960e-04 * zt + 0.378110e-02 ) * zs & & + ( ( - 0.302285e-13 * zh & & - 0.251520e-11 * zs & & + 0.512857e-12 * zt * zt ) * zh & & - 0.164759e-06 * zs & & +( 0.791325e-08 * zt - 0.933746e-06 ) * zt & & + 0.380374e-04 ) * zh zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta & - 0.301985e-05 ) * zt & & + 0.785567e-03 & & + ( 0.515032e-08 * zs & & + 0.788212e-08 * zt - 0.356603e-06 ) * zs & & +( ( 0.121551e-17 * zh & & - 0.602281e-15 * zs & & - 0.175379e-14 * zt + 0.176621e-12 ) * zh & & + 0.408195e-10 * zs & & + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt & & - 0.121555e-07 ) * zh !! tangent part zttl = 0.5 * ( ptem_tl(ji,jj,jk) + ptem_tl(ji,jj,jk-1) ) ! potential temperature at w-point zstl = 0.5 * ( psal_tl(ji,jj,jk) + psal_tl(ji,jj,jk-1) ) ! salinity anomaly at w-point zalbettl = ( ( ( -4.*0.255019e-07 * zt &! ratio alpha/beta & +3.*0.298357e-05 ) * zt & & -2.*0.203814e-03 ) * zt & & + 0.170907e-01 & & - 0.846960e-04 * zs & & - ( 0.933746e-06 & & - ( 2.*0.791325e-08 & & +2.*0.512857e-12 * zh ) * zt ) * zh ) * zttl & & + ( - 2.*0.678662e-05 * zs & & - 0.846960e-04 * zt & & + 0.378110e-02 & & + ( - 0.164759e-06 & & - 0.251520e-11 * zh ) * zh ) * zstl zbetatl = ( ( -3.*0.415613e-09 * zt & & +2.*0.555579e-07 ) * zt & & - 0.301985e-05 & & + 0.788212e-08 * zs & & + ( -2.*0.213127e-11 * zt & & - 0.175379e-14 * zh & & + 0.192867e-09 ) * zh ) * zttl & & + ( 2.*0.515032e-08 * zs & & + 0.788212e-08 * zt & & - 0.356603e-06 & & + ( - 0.602281e-15 * zh & & + 0.408195e-10 ) * zh ) * zstl pn2_tl(ji,jj,jk) = zgde3w * tmask(ji,jj,jk) * ( & & zbeta * ( zalbet & & * ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) & & + zalbettl & & * ( ptem (ji,jj,jk-1) - ptem (ji,jj,jk) ) & & - ( psal_tl(ji,jj,jk-1) - psal_tl(ji,jj,jk) ) ) & & + zbetatl * ( zalbet & & * ( ptem (ji,jj,jk-1) - ptem (ji,jj,jk) ) & & - ( psal (ji,jj,jk-1) - psal (ji,jj,jk) ) ) ) #if defined key_zdfddm zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) ) zdstl = ( psal_tl(ji,jj,jk-1) - psal_tl(ji,jj,jk) ) IF ( ABS( zds) <= 1.e-20 ) THEN zds = 1.e-20 rrau_tl(ji,jj,jk) = zalbettl * & & ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds & & + zalbet * & & ( ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) / zds ) ELSE rrau_tl(ji,jj,jk) = zalbettl * & & ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds & & + zalbet * & & ( ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) / zds & & - ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) * zdstl/ zds**2 ) ENDIF #endif END DO END DO END DO ! CASE ( 1 ) !== Linear formulation = F( temperature ) ==! DO jk = 2, jpkm1 pn2_tl(:,:,jk) = rn_alpha * ( ptem_tl(:,:,jk-1) - ptem_tl(:,:,jk) ) * grav / fse3w(:,:,jk) * tmask(:,:,jk) END DO ! CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! DO jk = 2, jpkm1 pn2_tl(:,:,jk) = ( rn_alpha * ( ptem_tl(:,:,jk-1) - ptem_tl(:,:,jk) ) & & - rn_beta * ( psal_tl(:,:,jk-1) - psal_tl(:,:,jk) ) ) & & * grav / fse3w(:,:,jk) * tmask(:,:,jk) END DO #if defined key_zdfddm DO jk = 2, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) ) zdstl = psal_tl(ji,jj,jk-1) - psal_tl(ji,jj,jk) IF ( ABS( zds) <= 1.e-20 ) THEN zds = 1.e-20 rrau_tl(ji,jj,jk) = ralpbet * & & ( ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) / zds ) ELSE rrau_tl(ji,jj,jk) = ralpbet * & & ( ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) / zds & & - ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) * zdstl / zds**2 ) ENDIF rrau(ji,jj,jk) = ralpbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds END DO END DO END DO #endif END SELECT END SUBROUTINE eos_bn2_tan SUBROUTINE eos_bn2_adj ( ptem, psal, ptem_ad, psal_ad, pn2_ad ) !!---------------------------------------------------------------------- !! *** ROUTINE eos_bn2_adj *** !! !! ** Purpose of the direct routine: Compute the local !! Brunt-Vaisala frequency at the time-step of the input arguments !! !! ** Method of the direct routine: !! * nn_eos = 0 : UNESCO sea water properties !! The brunt-vaisala frequency is computed using the polynomial !! polynomial expression of McDougall (1987): !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is !! computed and used in zdfddm module : !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) !! * nn_eos = 1 : linear equation of state (temperature only) !! N^2 = grav * rn_alpha * dk[ t ]/e3w !! * nn_eos = 2 : linear equation of state (temperature & salinity) !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w !! The use of potential density to compute N^2 introduces e r r o r !! in the sign of N^2 at great depths. We recommand the use of !! nn_eos = 0, except for academical studies. !! Macro-tasked on horizontal slab (jk-loop) !! N.B. N^2 is set to zero at the first level (JK=1) in inidtr !! and is never used at this level. !! !! ** Action : - pn2 : the brunt-vaisala frequency !! !! References : !! McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987. !! !! History: !! ! 08-07 (A. Vidard) First version !!---------------------------------------------------------------------- !! * Arguments REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & ptem, & ! potential temperature psal ! salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & ptem_ad, & ! adjoint potential temperature psal_ad ! adjoint salinity REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & pn2_ad ! adjoint Brunt-Vaisala frequency !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zgde3w, zt, zs, zh, & ! temporary scalars zalbet, zbeta ! " " REAL(wp) :: & ztad, zsad, & ! temporary scalars zalbetad, zbetaad ! " " #if defined key_zdfddm REAL(wp) :: zds, zdsad ! temporary scalars #endif ! pn2_tl : interior points only (2=< jk =< jpkm1 ) ! -------------------------- zalbetad = 0.0_wp zbetaad = 0.0_wp ztad = 0.0_wp zsad = 0.0_wp #if defined key_zdfddm zdsad = 0.0_wp #endif SELECT CASE ( nn_eos ) CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==! DO jk = jpkm1, 2, -1 DO jj = jpj, 1, -1 DO ji = jpi, 1, -1 zgde3w = grav / fse3w(ji,jj,jk) zt = 0.5 * ( ptem(ji,jj,jk) + ptem(ji,jj,jk-1) ) ! potential temperature at w-point zs = 0.5 * ( psal(ji,jj,jk) + psal(ji,jj,jk-1) ) - 35.0 ! salinity anomaly (s-35) at w-point zh = fsdepw(ji,jj,jk) ! depth in meters at w-point zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta & - 0.203814e-03 ) * zt & & + 0.170907e-01 ) * zt & & + 0.665157e-01 & & + ( - 0.678662e-05 * zs & & - 0.846960e-04 * zt + 0.378110e-02 ) * zs & & + ( ( - 0.302285e-13 * zh & & - 0.251520e-11 * zs & & + 0.512857e-12 * zt * zt ) * zh & & - 0.164759e-06 * zs & & +( 0.791325e-08 * zt - 0.933746e-06 ) * zt & & + 0.380374e-04 ) * zh zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta & - 0.301985e-05 ) * zt & & + 0.785567e-03 & & + ( 0.515032e-08 * zs & & + 0.788212e-08 * zt - 0.356603e-06 ) * zs & & +( ( 0.121551e-17 * zh & & - 0.602281e-15 * zs & & - 0.175379e-14 * zt + 0.176621e-12 ) * zh & & + 0.408195e-10 * zs & & + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt & & - 0.121555e-07 ) * zh #if defined key_zdfddm zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) ) IF ( ABS( zds) <= 1.e-20 ) THEN zds = 1.e-20 zdsad = 0.0_wp ELSE zdsad = rrau_ad(ji,jj,jk) * zalbet *( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds**2 ENDIF ptem_ad(ji,jj,jk-1) = ptem_ad(ji,jj,jk-1) + rrau_ad(ji,jj,jk) * zalbet / zds ptem_ad(ji,jj,jk ) = ptem_ad(ji,jj,jk ) - rrau_ad(ji,jj,jk) * zalbet / zds zalbetad = zalbetad + rrau_ad(ji,jj,jk) * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds rrau_ad(ji,jj,jk) = 0._wp psal_ad(ji,jj,jk-1) = psal_ad(ji,jj,jk-1) + zdsad psal_ad(ji,jj,jk ) = psal_ad(ji,jj,jk ) - zdsad zdsad = 0._wp #endif ptem_ad(ji,jj,jk-1) = ptem_ad(ji,jj,jk-1) + zalbet*zbeta*zgde3w*tmask(ji,jj,jk)*pn2_ad(ji,jj,jk) ptem_ad(ji,jj,jk) = ptem_ad(ji,jj,jk ) - zalbet*zbeta*zgde3w*tmask(ji,jj,jk)*pn2_ad(ji,jj,jk) zalbetad = zalbetad + zbeta*zgde3w*tmask(ji,jj,jk)*( ptem (ji,jj,jk-1) - ptem (ji,jj,jk) ) *pn2_ad(ji,jj,jk) psal_ad(ji,jj,jk-1) = psal_ad(ji,jj,jk-1) - zbeta*tmask(ji,jj,jk)*zgde3w*pn2_ad(ji,jj,jk) psal_ad(ji,jj,jk ) = psal_ad(ji,jj,jk ) + zbeta*tmask(ji,jj,jk)*zgde3w*pn2_ad(ji,jj,jk) zbetaad = zbetaad & & + zgde3w *tmask(ji,jj,jk)* ( zalbet * ( ptem (ji,jj,jk-1) - ptem (ji,jj,jk) ) & & - ( psal (ji,jj,jk-1) - psal (ji,jj,jk) ) )*pn2_ad(ji,jj,jk) pn2_ad(ji,jj,jk) = 0.0_wp ztad = ztad + ( ( -3.*0.415613e-09 * zt & & +2.*0.555579e-07 ) * zt & & - 0.301985e-05 & & + 0.788212e-08 * zs & & + ( -2.*0.213127e-11 * zt & & - 0.175379e-14 * zh & & + 0.192867e-09 ) * zh ) *zbetaad zsad = zsad + ( 2.*0.515032e-08 * zs & & + 0.788212e-08 * zt & & - 0.356603e-06 & & + ( - 0.602281e-15 * zh & & + 0.408195e-10 ) * zh ) * zbetaad zbetaad = 0.0_wp ztad = ztad + ( ( ( -4.*0.255019e-07 * zt &! ratio alpha/beta & +3.*0.298357e-05 ) * zt & & -2.*0.203814e-03 ) * zt & & + 0.170907e-01 & & - 0.846960e-04 * zs & & - ( 0.933746e-06 & & - ( 2.*0.791325e-08 & & +2.*0.512857e-12 * zh ) * zt ) * zh & & ) *zalbetad zsad = zsad + ( - 2.*0.678662e-05 * zs & & - 0.846960e-04 * zt & & + 0.378110e-02 & & + ( - 0.164759e-06 & & - 0.251520e-11 * zh ) * zh & & ) *zalbetad zalbetad = 0.0_wp psal_ad(ji,jj,jk) = psal_ad(ji,jj,jk) + 0.5 * zsad psal_ad(ji,jj,jk-1) = psal_ad(ji,jj,jk-1) + 0.5 * zsad zsad = 0.0_wp ptem_ad(ji,jj,jk) = ptem_ad(ji,jj,jk) + 0.5 * ztad ptem_ad(ji,jj,jk-1) = ptem_ad(ji,jj,jk-1) + 0.5 * ztad ztad = 0.0_wp END DO END DO END DO ! CASE ( 1 ) !== Linear formulation = F( temperature ) ==! DO jk = jpkm1, 2, -1 ptem_ad(:,:,jk-1) = ptem_ad(:,:,jk-1) + rn_alpha * pn2_ad(:,:,jk) & & * grav / fse3w(:,:,jk) * tmask(:,:,jk) ptem_ad(:,:,jk ) = ptem_ad(:,:,jk ) - rn_alpha * pn2_ad(:,:,jk) & & * grav / fse3w(:,:,jk) * tmask(:,:,jk) pn2_ad(:,:,jk) = 0.0_wp END DO ! CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! #if defined key_zdfddm DO jk = jpkm1, 2, -1 DO jj = jpj, 1, -1 DO ji = jpi, 1, -1 zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) ) IF ( ABS( zds) <= 1.e-20 ) THEN zds = 1.e-20 zdsad = 0.0_wp ELSE zdsad = zdsad - rrau_ad(ji,jj,jk) * ralpbet & & * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds**2 ENDIF rrau(ji,jj,jk) = ralpbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds ptem_ad(ji,jj,jk-1) = ptem_ad(ji,jj,jk-1) & & + rrau_ad(ji,jj,jk) * ralpbet / zds ptem_ad(ji,jj,jk ) = ptem_ad(ji,jj,jk ) & & - rrau_ad(ji,jj,jk) * ralpbet / zds rrau_ad(ji,jj,jk) = 0._wp psal_ad(ji,jj,jk-1) = psal_ad(ji,jj,jk-1) + zdsad psal_ad(ji,jj,jk ) = psal_ad(ji,jj,jk ) - zdsad zdsad = 0._wp END DO END DO END DO #endif DO jk = jpkm1, 2, -1 ptem_ad(:,:,jk-1) = ptem_ad(:,:,jk-1) + rn_alpha * pn2_ad(:,:,jk) & & * grav / fse3w(:,:,jk) * tmask(:,:,jk) ptem_ad(:,:,jk ) = ptem_ad(:,:,jk ) - rn_alpha * pn2_ad(:,:,jk) & & * grav / fse3w(:,:,jk) * tmask(:,:,jk) psal_ad(:,:,jk-1) = psal_ad(:,:,jk-1) - rn_beta * pn2_ad(:,:,jk) & & * grav / fse3w(:,:,jk) * tmask(:,:,jk) psal_ad(:,:,jk ) = psal_ad(:,:,jk ) + rn_beta * pn2_ad(:,:,jk) & & * grav / fse3w(:,:,jk) * tmask(:,:,jk) pn2_ad(:,:,jk) = 0.0_wp END DO END SELECT END SUBROUTINE eos_bn2_adj #if defined key_tam SUBROUTINE eos_insitu_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_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 !! !! !! History : !! ! 08-07 (A. Vidard) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & ztem, & ! potential temperature zsal ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zt_adout, & ! potential temperature & zs_adout ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrd_adin ! anomaly density REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zt_tlin, & ! potential temperature & zs_tlin ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & znt, & ! potential temperature & zns ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrd_tlout ! anomaly density REAL(KIND=wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp2, & ! scalar product involving the adjoint routine & zsp2_1, & ! scalar product involving the adjoint routine & zsp2_2 ! scalar product involving the adjoint routine INTEGER, DIMENSION(jpi,jpj) :: & & iseed_2d ! 2D seed for the random number generator INTEGER :: & & ji, & & jj, & & jk, & & jn, & & jeos CHARACTER(LEN=14) :: cl_name ALLOCATE( & & ztem( jpi, jpj, jpk ), & & zsal( jpi, jpj, jpk ), & & znt( jpi, jpj, jpk ), & & zns( jpi, jpj, jpk ), & & zt_adout( jpi, jpj, jpk ), & & zs_adout( jpi, jpj, jpk ), & & zrd_adin( jpi, jpj, jpk ), & & zs_tlin( jpi, jpj, jpk ), & & zt_tlin( jpi, jpj, jpk ), & & zrd_tlout(jpi, jpj, jpk ) ) ! Initialize the reference state ztem(:,:,:) = tn(:,:,:) zsal(:,:,:) = sn(:,:,:) ! store initial nn_eos jeos = nn_eos DO jn = 0, 2 nn_eos = jn !============================================================= ! 1) dx = ( T ) and dy = ( T ) !============================================================= !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zt_tlin(:,:,:) = 0.0_wp zs_tlin(:,:,:) = 0.0_wp zrd_tlout(:,:,:) = 0.0_wp zt_adout(:,:,:) = 0.0_wp zs_adout(:,:,:) = 0.0_wp zrd_adin(:,:,:) = 0.0_wp !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- 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, znt, 'T', 0.0_wp, stdt ) DO jj = 1, jpj DO ji = 1, jpi iseed_2d(ji,jj) = - ( 395703 + & & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) END DO END DO CALL grid_random( iseed_2d, zns, 'T', 0.0_wp, stds ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zt_tlin(ji,jj,jk) = znt(ji,jj,jk) zs_tlin(ji,jj,jk) = zns(ji,jj,jk) END DO END DO END DO CALL eos_insitu_tan(ztem, zsal, zt_tlin, zs_tlin, zrd_tlout) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zrd_adin(ji,jj,jk) = zrd_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( zrd_tlout, zrd_adin ) !-------------------------------------------------------------------- ! Call the adjoint routine: dx^* = L^T dy^* !-------------------------------------------------------------------- CALL eos_insitu_adj(ztem, zsal, zt_adout, zs_adout, zrd_adin) zsp2_1 = DOT_PRODUCT( zt_tlin, zt_adout ) zsp2_2 = DOT_PRODUCT( zs_tlin, zs_adout ) zsp2 = zsp2_1 + zsp2_2 ! Compare the scalar products ! Compare the scalar products ! 14 char:'12345678901234' SELECT CASE( jn ) CASE (0) ; cl_name = 'eos_adj ins T1' CASE (1) ; cl_name = 'eos_adj ins T2' CASE (2) ; cl_name = 'eos_adj ins T3' END SELECT CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) ENDDO ! restore initial nn_eos nn_eos = jeos ! Deallocate memory DEALLOCATE( & & ztem, & & zsal, & & zt_adout, & & zs_adout, & & zrd_adin, & & zt_tlin, & & zs_tlin, & & zrd_tlout, & & znt, & & zns & & ) END SUBROUTINE eos_insitu_adj_tst SUBROUTINE eos_insitu_pot_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_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-07 (A. Vidard) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & ztem, & ! potential temperature zsal ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zt_adout, & ! potential temperature & zs_adout ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrd_adin ! anomaly density REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrhop_adin ! volume mass REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zt_tlin, & ! potential temperature & zs_tlin ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & znt, & ! potential temperature & zns ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrd_tlout ! anomaly density REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrhop_tlout ! volume mass REAL(KIND=wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp2, & ! scalar product involving the adjoint routine & zsp2_1, & ! scalar product involving the adjoint routine & zsp2_2 ! scalar product involving the adjoint routine INTEGER, DIMENSION(jpi,jpj) :: & & iseed_2d ! 2D seed for the random number generator INTEGER :: & & ji, & & jj, & & jk, & & jn, & & jeos CHARACTER(LEN=14) :: cl_name ! Allocate memory ALLOCATE( & & ztem( jpi, jpj, jpk ), & & zsal( jpi, jpj, jpk ), & & zt_adout( jpi, jpj, jpk ), & & zs_adout( jpi, jpj, jpk ), & & zrhop_adin( jpi, jpj, jpk ), & & zrd_adin( jpi, jpj, jpk ), & & zs_tlin( jpi, jpj, jpk ), & & zt_tlin( jpi, jpj, jpk ), & & zns( jpi, jpj, jpk ), & & znt( jpi, jpj, jpk ), & & zrd_tlout(jpi, jpj, jpk ), & & zrhop_tlout(jpi, jpj, jpk ) ) ! Initialize random field standard deviationsthe reference state ztem = tn zsal = sn ! store initial nn_eos jeos = nn_eos DO jn = 0, 2 nn_eos = jn !============================================= ! testing of eos_insitu_pot !============================================= !============================================================= ! 1) dx = ( T ) and dy = ( T ) !============================================================= !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zt_tlin(:,:,:) = 0.0_wp zs_tlin(:,:,:) = 0.0_wp zrd_tlout(:,:,:) = 0.0_wp zrhop_tlout(:,:,:) = 0.0_wp zt_adout(:,:,:) = 0.0_wp zs_adout(:,:,:) = 0.0_wp zrhop_adin(:,:,:) = 0.0_wp zrd_adin(:,:,:) = 0.0_wp !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- 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, znt, 'T', 0.0_wp, stdt ) 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, zns, 'T', 0.0_wp, stds ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zt_tlin(ji,jj,jk) = znt(ji,jj,jk) zs_tlin(ji,jj,jk) = zns(ji,jj,jk) END DO END DO END DO !-------------------------------------------------------------------- ! Call the tangent routine: dy = L dx !-------------------------------------------------------------------- call eos_insitu_pot_tan ( ztem, zsal, zt_tlin, zs_tlin, zrd_tlout, zrhop_tlout ) !-------------------------------------------------------------------- ! Initialize the adjoint variables: dy^* = W dy !-------------------------------------------------------------------- DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) & & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)& & * tmask(ji,jj,jk) zrhop_adin(ji,jj,jk) = zrhop_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( zrd_tlout , zrd_adin ) & & + DOT_PRODUCT( zrhop_tlout, zrhop_adin ) !-------------------------------------------------------------------- ! Call the adjoint routine: dx^* = L^T dy^* !-------------------------------------------------------------------- CALL eos_insitu_pot_adj( ztem, zsal, zt_adout, zs_adout, zrd_adin, zrhop_adin ) !-------------------------------------------------------------------- ! Compute the scalar product: dx^T L^T W dy !-------------------------------------------------------------------- zsp2_1 = DOT_PRODUCT( zt_tlin, zt_adout ) zsp2_2 = DOT_PRODUCT( zs_tlin, zs_adout ) zsp2 = zsp2_1 + zsp2_2 ! Compare the scalar products ! 14 char:'12345678901234' SELECT CASE( jn ) CASE (0) ; cl_name = 'eos_adj pot T1' CASE (1) ; cl_name = 'eos_adj pot T2' CASE (2) ; cl_name = 'eos_adj pot T3' END SELECT CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) ENDDO ! restore initial nn_eos nn_eos = jeos ! Deallocate memory DEALLOCATE( & & ztem, & & zsal, & & zt_adout, & & zs_adout, & & zrd_adin, & & zrhop_adin, & & zt_tlin, & & zs_tlin, & & zrd_tlout, & & zrhop_tlout,& & zns, znt ) END SUBROUTINE eos_insitu_pot_adj_tst SUBROUTINE eos_insitu_2d_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_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-07 (A. Vidard) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & zdep ! depth REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & ztem, & ! potential temperature zsal ! salinity REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & zt_adout, & ! potential temperature & zs_adout ! salinity REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & zrd_adin ! anomaly density REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & zt_tlin, & ! potential temperature & zs_tlin ! salinity REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & znt, & ! potential temperature & zns ! salinity REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & zrd_tlout ! anomaly density REAL(KIND=wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp2, & ! scalar product involving the adjoint routine & zsp2_1, & ! scalar product involving the adjoint routine & zsp2_2 ! scalar product involving the adjoint routine INTEGER, DIMENSION(jpi,jpj) :: & & iseed_2d ! 2D seed for the random number generator INTEGER :: & & ji, & & jj, & & jn, & & jeos CHARACTER(LEN=14) :: cl_name ! Allocate memory ALLOCATE( & & zdep( jpi, jpj ), & & ztem( jpi, jpj ), & & zsal( jpi, jpj ), & & znt( jpi, jpj ), & & zns( jpi, jpj ), & & zt_adout( jpi, jpj ), & & zs_adout( jpi, jpj ), & & zrd_adin( jpi, jpj ), & & zs_tlin( jpi, jpj ), & & zt_tlin( jpi, jpj ), & & zrd_tlout(jpi, jpj ) ) ! Initialize the reference state ztem(:,:) = tn(:,:,2) zsal(:,:) = sn(:,:,2) zdep(:,:) = fsdept(:,:,2) ! store initial nn_eos jeos = nn_eos DO jn = 0, 2 nn_eos = jn !============================================================= ! 1) dx = ( T ) and dy = ( T ) !============================================================= !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zt_tlin(:,:) = 0.0_wp zs_tlin(:,:) = 0.0_wp zrd_tlout(:,:) = 0.0_wp zt_adout(:,:) = 0.0_wp zs_adout(:,:) = 0.0_wp zrd_adin(:,:) = 0.0_wp !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- 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, znt, 'T', 0.0_wp, stdt ) 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, zns, 'T', 0.0_wp, stds ) DO jj = nldj, nlej DO ji = nldi, nlei zt_tlin(ji,jj) = znt(ji,jj) zs_tlin(ji,jj) = zns(ji,jj) END DO END DO CALL eos_insitu_2d_tan(ztem, zsal, zdep, zt_tlin, zs_tlin, zrd_tlout) DO jj = nldj, nlej DO ji = nldi, nlei zrd_adin(ji,jj) = zrd_tlout(ji,jj) & & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,2)& & * tmask(ji,jj,2) END DO END DO !-------------------------------------------------------------------- ! Compute the scalar product: ( L dx )^T W dy !-------------------------------------------------------------------- zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin ) !-------------------------------------------------------------------- ! Call the adjoint routine: dx^* = L^T dy^* !-------------------------------------------------------------------- CALL eos_insitu_2d_adj(ztem, zsal, zdep, zt_adout, zs_adout, zrd_adin) zsp2_1 = DOT_PRODUCT( zt_tlin, zt_adout ) zsp2_2 = DOT_PRODUCT( zs_tlin, zs_adout ) zsp2 = zsp2_1 + zsp2_2 ! Compare the scalar products ! 14 char:'12345678901234' SELECT CASE( jn ) CASE (0) ; cl_name = 'eos_adj 2d T1' CASE (1) ; cl_name = 'eos_adj 2d T2' CASE (2) ; cl_name = 'eos_adj 2d T3' END SELECT CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) ENDDO ! restore initial nn_eos nn_eos = jeos ! Deallocate memory DEALLOCATE( & & zdep, & & ztem, & & zsal, & & zt_adout, & & zs_adout, & & zrd_adin, & & zt_tlin, & & zs_tlin, & & zrd_tlout, & & zns, znt ) END SUBROUTINE eos_insitu_2d_adj_tst SUBROUTINE eos_pot_1pt_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_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-07 (A. Vidard) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit !! * Local declarations REAL(KIND=wp) :: & & ztem, & & zsal, & & zt_tlin, & & zs_tlin, & & zrhop_tlout, & & zt_adout, & & zs_adout, & & zrhop_adin REAL(KIND=wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp2 ! scalar product involving the adjoint routine CHARACTER(LEN=14) :: cl_name INTEGER :: & & jn, & & jeos ! Initialize the reference state ztem = 23.7 zsal = 30.1 ! store initial nn_eos jeos = nn_eos DO jn = 0, 2 nn_eos = jn !================================================================== ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and ! dy = ( hdivb_tl, hdivn_tl ) !================================================================== !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zt_tlin = 1.12_wp zs_tlin = 0.123_wp zrhop_tlout = 0.0_wp zt_adout = 0.0_wp zs_adout = 0.0_wp zrhop_adin = 0.0_wp CALL eos_pot_1pt_tan( ztem, zsal, zt_tlin, zs_tlin, zrhop_tlout ) !-------------------------------------------------------------------- ! Initialize the adjoint variables: dy^* = W dy !-------------------------------------------------------------------- zrhop_adin = zrhop_tlout * e1t(1,1) * e2t(1,1) * fse3t(1,1,1) !-------------------------------------------------------------------- ! Compute the scalar product: ( L dx )^T W dy !-------------------------------------------------------------------- zsp1 = zrhop_adin * zrhop_tlout CALL eos_pot_1pt_adj( ztem, zsal, zt_adout, zs_adout, zrhop_adin ) zsp2 = zt_tlin * zt_adout + zs_tlin * zs_adout ! 14 char:'12345678901234' SELECT CASE( jn ) CASE (0) ; cl_name = 'eos_adj 1pt T1' CASE (1) ; cl_name = 'eos_adj 1pt T2' CASE (2) ; cl_name = 'eos_adj 1pt T3' END SELECT CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) ENDDO ! restore initial nn_eos nn_eos = jeos END SUBROUTINE eos_pot_1pt_adj_tst SUBROUTINE eos_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_adj_tst *** !! !! ** Purpose : Test the adjoint routine. !! !! History : !! ! 08-07 (A. Vidard) !!----------------------------------------------------------------------- !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit CALL eos_insitu_adj_tst( kumadt ) CALL eos_insitu_pot_adj_tst( kumadt ) CALL eos_insitu_2d_adj_tst( kumadt ) CALL eos_pot_1pt_adj_tst( kumadt ) END SUBROUTINE eos_adj_tst SUBROUTINE bn2_adj_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE bn2_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-07 (A. Vidard) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & ztem, & ! potential temperature zsal ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zt_adout, & ! potential temperature & zs_adout ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrd_adin, & ! potential density (surface referenced) & zrd_adout ! potential density (surface referenced) REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zt_tlin, & ! potential temperature & zs_tlin, & ! salinity & zt_tlout, & ! potential temperature & zs_tlout ! salinity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrd_tlout ! potential density (surface referenced) REAL(KIND=wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp2, & ! scalar product involving the adjoint routine & zsp2_1, & ! scalar product involving the adjoint routine & zsp2_2 ! scalar product involving the adjoint routine REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & znt, & ! potential temperature & zns ! salinity INTEGER, DIMENSION(jpi,jpj) :: & & iseed_2d ! 2D seed for the random number generator INTEGER :: & & iseed, & & ji, & & jj, & & jk, & & jn, & & jeos CHARACTER(LEN=14) :: cl_name ! Allocate memory ALLOCATE( & & ztem( jpi, jpj, jpk ), & & zsal( jpi, jpj, jpk ), & & zt_adout( jpi, jpj, jpk ), & & zs_adout( jpi, jpj, jpk ), & & zrd_adin( jpi, jpj, jpk ), & & zrd_adout(jpi, jpj, jpk ), & & zs_tlin( jpi, jpj, jpk ), & & zt_tlin( jpi, jpj, jpk ), & & zns( jpi, jpj, jpk ), & & znt( jpi, jpj, jpk ), & & zt_tlout( jpi, jpj, jpk ), & & zs_tlout( jpi, jpj, jpk ), & & zrd_tlout(jpi, jpj, jpk ) ) ! Initialize random field standard deviationsthe reference state ztem = tn zsal = sn ! store initial nn_eos jeos = nn_eos DO jn = 0, 2 nn_eos = jn !============================================================= ! 1) dx = ( T ) and dy = ( T ) !============================================================= !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zt_tlin(:,:,:) = 0.0_wp zs_tlin(:,:,:) = 0.0_wp zt_tlout(:,:,:) = 0.0_wp zs_tlout(:,:,:) = 0.0_wp zrd_tlout(:,:,:) = 0.0_wp zt_adout(:,:,:) = 0.0_wp zs_adout(:,:,:) = 0.0_wp zrd_adin(:,:,:) = 0.0_wp zrd_adout(:,:,:) = 0.0_wp !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- 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, znt, 'T', 0.0_wp, stdt ) 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, zns, 'T', 0.0_wp, stds ) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zt_tlin(ji,jj,jk) = znt(ji,jj,jk) zs_tlin(ji,jj,jk) = zns(ji,jj,jk) END DO END DO END DO !-------------------------------------------------------------------- ! Call the tangent routine: dy = L dx !-------------------------------------------------------------------- zt_tlout(:,:,:) = zt_tlin zs_tlout(:,:,:) = zs_tlin CALL eos_bn2_tan( ztem, zsal, zt_tlout, zs_tlout, zrd_tlout ) !-------------------------------------------------------------------- ! Initialize the adjoint variables: dy^* = W dy !-------------------------------------------------------------------- DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zrd_adin(ji,jj,jk) = zrd_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( zrd_tlout, zrd_adin ) !-------------------------------------------------------------------- ! Call the adjoint routine: dx^* = L^T dy^* !-------------------------------------------------------------------- zrd_adout(:,:,:) = zrd_adin(:,:,:) CALL eos_bn2_adj( ztem, zsal, zt_adout, zs_adout, zrd_adout ) !-------------------------------------------------------------------- ! Compute the scalar product: dx^T L^T W dy !-------------------------------------------------------------------- zsp2_1 = DOT_PRODUCT( zt_tlin, zt_adout ) zsp2_2 = DOT_PRODUCT( zs_tlin, zs_adout ) zsp2 = zsp2_1 + zsp2_2 ! Compare the scalar products ! 14 char:'12345678901234' SELECT CASE( jn ) CASE (0) ; cl_name = 'bn2_adj T1' CASE (1) ; cl_name = 'bn2_adj T2' CASE (2) ; cl_name = 'bn2_adj T3' END SELECT CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) ENDDO ! restore initial nn_eos nn_eos = jeos ! Deallocate memory DEALLOCATE( & & ztem, & & zsal, & & zt_adout, & & zs_adout, & & zrd_adin, & & zrd_adout, & & zt_tlin, & & zs_tlin, & & zt_tlout, & & zs_tlout, & & zrd_tlout, & & zns, znt ) END SUBROUTINE bn2_adj_tst #if defined key_tst_tlm SUBROUTINE eos_insitu_tlm_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_insitu_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 eosbn2, ONLY: & ! horizontal & vertical advective trend & eos 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 & tn, sn, rhd, rhop USE in_out_manager, ONLY: & ! I/O manager & nitend, & & nit000 USE tamctl, ONLY: & ! Control parameters & numtan, numtan_sc !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrd_out, & ! & zt_tlin , & ! & zs_tlin , & & zrd_tl , & & zrd_wop, & & z3r REAL(KIND=wp) :: & & zsp1, & & zsp2, & & zsp3, & & zzsp, & & gamma, & & zgsp1, & & zgsp2, & & zgsp3, & & zgsp4, & & zgsp5, & & zgsp6, & & zgsp7 INTEGER :: & & ji, & & jj, & & jk CHARACTER(LEN=14) :: cl_name CHARACTER (LEN=128) :: file_out_sc, file_wop, file_out, file_xdx CHARACTER (LEN=90) :: FMT REAL(KIND=wp), DIMENSION(100):: & & zscrd, zscerrrd INTEGER, DIMENSION(100):: & & iiposrd, ijposrd, ikposrd INTEGER:: & & ii, numsctlm, & & numtlm, & & isamp=40,jsamp=40, ksamp=10 REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: & & zerrrd ALLOCATE( & & zrd_out( jpi, jpj, jpk ), & & zrd_tl( jpi, jpj, jpk ), & & zs_tlin( jpi, jpj, jpk ), & & zt_tlin( jpi, jpj, jpk ), & & zrd_wop( jpi, jpj, jpk ), & & z3r ( jpi, jpj, jpk ) ) !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zt_tlin( :,:,:) = 0.0_wp zs_tlin( :,:,:) = 0.0_wp zrd_out( :,:,:) = 0.0_wp zrd_wop( :,:,:) = 0.0_wp zscerrrd(:) = 0.0_wp zscrd(:) = 0.0_wp IF ( tlm_bch == 2 ) zrd_tl ( :,:,:) = 0.0_wp !-------------------------------------------------------------------- ! Output filename Xn=F(X0) !-------------------------------------------------------------------- ! CALL tlm_namrd gamma = h_ratio file_wop='trj_wop_eos_insitu' file_xdx='trj_xdx_eos_insitu' !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN CALL grid_rd_sd( 596035, z3r, 'T', 0.0_wp, stdt) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zt_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 371836, z3r, 'S', 0.0_wp, stds) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zs_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO ENDIF !-------------------------------------------------------------------- ! Complete Init for Direct !------------------------------------------------------------------- IF ( tlm_bch /= 2 ) 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 zt_tlin(:,:,:) = gamma * zt_tlin(:,:,:) tn(:,:,:) = tn(:,:,:) + zt_tlin(:,:,:) zs_tlin(:,:,:) = gamma * zs_tlin(:,:,:) sn(:,:,:) = sn(:,:,:) + zs_tlin(:,:,:) ENDIF !-------------------------------------------------------------------- ! Compute the direct model F(X0,t=n) = Xn !-------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL eos(tn, sn, zrd_out) rhd(:,:,:)= zrd_out(:,:,:) 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: dy^* = W dy !-------------------------------------------------------------------- CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) !----------------------------------------------------------------------- ! Initialization of the dynamics and tracer fields for the tangent !----------------------------------------------------------------------- CALL eos_insitu_tan(tn, sn, zt_tlin, zs_tlin, zrd_tl) !-------------------------------------------------------------------- ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) !-------------------------------------------------------------------- zsp2 = DOT_PRODUCT( zrd_tl, zrd_tl ) !-------------------------------------------------------------------- ! Storing data !-------------------------------------------------------------------- CALL trj_rd_spl(file_wop) zrd_wop (:,:,:) = rhd (:,:,:) CALL trj_rd_spl(file_xdx) zrd_out (:,:,:) = rhd (:,:,:) !-------------------------------------------------------------------- ! 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 zrd_out (ji,jj,jk) = zrd_out (ji,jj,jk) - zrd_wop (ji,jj,jk) zrd_wop (ji,jj,jk) = zrd_out (ji,jj,jk) - zrd_tl (ji,jj,jk) IF ( zrd_tl(ji,jj,jk) .NE. 0.0_wp ) & & zerrrd(ji,jj,jk) = zrd_out(ji,jj,jk)/zrd_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 iiposrd(ii) = ji ijposrd(ii) = jj ikposrd(ii) = jk IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN zscrd (ii) = zrd_wop(ji,jj,jk) zscerrrd (ii) = ( zerrrd(ji,jj,jk) - 1.0_wp ) / gamma ENDIF ENDIF END DO END DO END DO zsp1 = DOT_PRODUCT( zrd_out, zrd_out ) zsp3 = DOT_PRODUCT( zrd_wop, zrd_wop ) !-------------------------------------------------------------------- ! Print the linearization error En - norme 2 !-------------------------------------------------------------------- ! 14 char:'12345678901234' cl_name = 'eos_insitu:En ' zzsp = SQRT(zsp3) zgsp5 = zzsp CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Compute TLM norm2 !-------------------------------------------------------------------- zzsp = SQRT(zsp2) zgsp4 = zzsp cl_name = 'eos_insitu:Ln2' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Print the linearization error Nn - norme 2 !-------------------------------------------------------------------- zzsp = SQRT(zsp1) cl_name = 'eosins: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) 'eosinsit', 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 = 'eosinsit' IF(lwp) THEN DO ii=1, 100, 1 IF ( zscrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrd ', & & cur_loop, h_ratio, iiposrd(ii), ijposrd(ii), ikposrd(ii), zscrd(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrrd', & & cur_loop, h_ratio, iiposrd(ii), ijposrd(ii), ikposrd(ii), zscerrrd(ii) ENDDO ! write separator WRITE(numtan_sc,"(A4)") '====' ENDIF ENDIF DEALLOCATE( & & zrd_out, zrd_tl, zrd_wop, & & zt_tlin, zs_tlin, z3r & & ) END SUBROUTINE eos_insitu_tlm_tst SUBROUTINE eos_insitu_pot_tlm_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_insitu_pot_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 eosbn2, ONLY: & ! horizontal & vertical advective trend & eos 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 & tn, sn, rhd, rhop USE in_out_manager, ONLY: & ! I/O manager & nitend, & & nit000 USE tamctl, ONLY: & ! Control parameters & numtan, numtan_sc !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & zrd_out, & ! & zrh_out, & & zt_tlin , & ! & zs_tlin , & & zrh_tl , & & zrd_tl , & & zrd_wop , & & zrh_wop , & & z3r REAL(KIND=wp) :: & & zsp1, & & zsp2, & & zsp3, & & zsp1_Rd, Zsp1_Rh, & & zsp2_Rd, zsp2_Rh, & & zsp3_Rd, zsp3_Rh, & & zzsp, & & gamma, & & zgsp1, & & zgsp2, & & zgsp3, & & zgsp4, & & zgsp5, & & zgsp6, & & zgsp7 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):: & & zscrd, zscerrrd, & & zscrh, zscerrrh INTEGER, DIMENSION(100):: & & iiposrd, ijposrd, ikposrd, & & iiposrh, ijposrh, ikposrh INTEGER:: & & ii, numsctlm, & & isamp=40,jsamp=40, ksamp=10 REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: & & zerrrd, zerrrh ALLOCATE( & & zrd_out( jpi, jpj, jpk ), & & zrh_out( jpi, jpj, jpk ), & & zrd_tl( jpi, jpj, jpk ), & & zrh_tl( jpi, jpj, jpk ), & & zs_tlin( jpi, jpj, jpk ), & & zt_tlin( jpi, jpj, jpk ), & & zrd_wop( jpi, jpj, jpk ), & & zrh_wop( jpi, jpj, jpk ), & & z3r ( jpi, jpj, jpk ) ) !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- zt_tlin( :,:,:) = 0.0_wp zs_tlin( :,:,:) = 0.0_wp zrd_out( :,:,:) = 0.0_wp zrh_out( :,:,:) = 0.0_wp zrd_wop( :,:,:) = 0.0_wp zrh_wop( :,:,:) = 0.0_wp zscerrrd( :) = 0.0_wp zscerrrh( :) = 0.0_wp zscrd(:) = 0.0_wp zscrh(:) = 0.0_wp IF ( tlm_bch == 2 ) THEN zrd_tl ( :,:,:) = 0.0_wp zrh_tl ( :,:,:) = 0.0_wp ENDIF !-------------------------------------------------------------------- ! Output filename Xn=F(X0) !-------------------------------------------------------------------- !! CALL tlm_namrd gamma = h_ratio file_wop='trj_wop_eos_pot' file_xdx='trj_xdx_eos_pot' !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN CALL grid_rd_sd( 596035, z3r, 'T', 0.0_wp, stdt) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zt_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 371836, z3r, 'S', 0.0_wp, stds) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zs_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO ENDIF !-------------------------------------------------------------------- ! Complete Init for Direct !------------------------------------------------------------------- IF ( tlm_bch /= 2 ) 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 zt_tlin(:,:,:) = gamma * zt_tlin(:,:,:) tn(:,:,:) = tn(:,:,:) + zt_tlin(:,:,:) zs_tlin(:,:,:) = gamma * zs_tlin(:,:,:) sn(:,:,:) = sn(:,:,:) + zs_tlin(:,:,:) ENDIF !-------------------------------------------------------------------- ! Compute the direct model F(X0,t=n) = Xn !-------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL eos(tn, sn, zrd_out, zrh_out) rhd (:,:,:) = zrd_out(:,:,:) rhop(:,:,:) = zrh_out(:,:,:) 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: dy^* = W dy !-------------------------------------------------------------------- CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) !----------------------------------------------------------------------- ! Initialization of the dynamics and tracer fields for the tangent !----------------------------------------------------------------------- CALL eos_insitu_pot_tan(tn, sn, zt_tlin, zs_tlin, zrd_tl, zrh_tl) !-------------------------------------------------------------------- ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) !-------------------------------------------------------------------- zsp2_Rd = DOT_PRODUCT( zrd_tl, zrd_tl ) zsp2_Rh = DOT_PRODUCT( zrh_tl, zrh_tl ) zsp2 = zsp2_Rd + zsp2_Rh !-------------------------------------------------------------------- ! Storing data !-------------------------------------------------------------------- CALL trj_rd_spl(file_wop) zrd_wop (:,:,:) = rhd (:,:,:) zrh_wop (:,:,:) = rhop (:,:,:) CALL trj_rd_spl(file_xdx) zrd_out (:,:,:) = rhd (:,:,:) zrh_out (:,:,:) = rhop (:,:,:) !-------------------------------------------------------------------- ! 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 zrd_out (ji,jj,jk) = zrd_out (ji,jj,jk) - zrd_wop (ji,jj,jk) zrd_wop (ji,jj,jk) = zrd_out (ji,jj,jk) - zrd_tl (ji,jj,jk) IF ( zrd_tl(ji,jj,jk) .NE. 0.0_wp ) & & zerrrd(ji,jj,jk) = zrd_out(ji,jj,jk)/zrd_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 iiposrd(ii) = ji ijposrd(ii) = jj ikposrd(ii) = jk IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN zscrd (ii) = zrd_wop(ji,jj,jk) zscerrrd (ii) = ( zerrrd( 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 zrh_out (ji,jj,jk) = zrh_out (ji,jj,jk) - zrh_wop (ji,jj,jk) zrh_wop (ji,jj,jk) = zrh_out (ji,jj,jk) - zrh_tl (ji,jj,jk) IF ( zrh_tl(ji,jj,jk) .NE. 0.0_wp ) & & zerrrh(ji,jj,jk) = zrh_out(ji,jj,jk)/zrh_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 iiposrh(ii) = ji ijposrh(ii) = jj ikposrh(ii) = jk IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN zscrh (ii) = zrh_wop(ji,jj,jk) zscerrrh (ii) = ( zerrrh( ji,jj,jk) - 1.0_wp ) /gamma ENDIF ENDIF END DO END DO END DO zsp1_Rd = DOT_PRODUCT( zrd_out, zrd_out ) zsp1_Rh = DOT_PRODUCT( zrh_out, zrh_out ) zsp1 = zsp1_Rd + zsp1_Rh zsp3_Rd = DOT_PRODUCT( zrd_wop, zrd_wop ) zsp3_Rh = DOT_PRODUCT( zrh_wop, zrh_wop ) zsp3 = zsp3_Rd + zsp3_Rh !-------------------------------------------------------------------- ! Print the linearization error En - norme 2 !-------------------------------------------------------------------- ! 14 char:'12345678901234' cl_name = 'eos_pot:En ' zzsp = SQRT(zsp3) zgsp5 = zzsp CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Compute TLM norm2 !-------------------------------------------------------------------- zzsp = SQRT(zsp2) zgsp4 = zzsp cl_name = 'eos_pot:Ln2' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Print the linearization error Nn - norme 2 !-------------------------------------------------------------------- zzsp = SQRT(zsp1) cl_name = 'eospot: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) 'eosinsp ', 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 = 'eosinsp ' IF(lwp) THEN DO ii=1, 100, 1 IF ( zscrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrd ', & & cur_loop, h_ratio, iiposrd(ii), ijposrd(ii), ikposrd(ii), zscrd(ii) ENDDO DO ii=1, 100, 1 IF ( zscrh(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrh ', & & cur_loop, h_ratio, iiposrh(ii), ijposrh(ii), ikposrh(ii), zscrh(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrrd', & & cur_loop, h_ratio, iiposrd(ii), ijposrd(ii), ikposrd(ii), zscerrrd(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrrh(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrrh', & & cur_loop, h_ratio, iiposrh(ii), ijposrh(ii), ikposrh(ii), zscerrrh(ii) ENDDO ! write separator WRITE(numtan_sc,"(A4)") '====' ENDIF ENDIF DEALLOCATE( & & zrd_out, zrd_tl, zrd_wop, & & zrh_out, zrh_tl, zrh_wop, & & zt_tlin, zs_tlin, z3r & & ) END SUBROUTINE eos_insitu_pot_tlm_tst SUBROUTINE eos_insitu_2d_tlm_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_insitu_2d_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 eosbn2, ONLY: & ! horizontal & vertical advective trend & eos 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 & tn, sn, rhd, rhop USE in_out_manager, ONLY: & ! I/O manager & nitend, & & nit000 USE tamctl, ONLY: & ! Control parameters & numtan, numtan_sc !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & zdep, & ! depth & ztem, & ! potential temperature & zsal ! salinity REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & & zrd_out, & ! & zt_tlin , & ! & zs_tlin , & & zrd_tl , & & zrd_wop , & & z2r REAL(KIND=wp) :: & & zsp1, & & zsp2, & & zsp3, & & zzsp, & & gamma, & & zgsp1, & & zgsp2, & & zgsp3, & & zgsp4, & & zgsp5, & & zgsp6, & & zgsp7 INTEGER :: & & ji, & & jj CHARACTER(LEN=14) :: cl_name CHARACTER (LEN=128) :: file_out, file_wop, file_xdx CHARACTER (LEN=90) :: FMT REAL(KIND=wp), DIMENSION(100):: & & zscrd, zscerrrd INTEGER, DIMENSION(100):: & & iiposrd, ijposrd INTEGER:: & & ii, numsctlm, & & isamp=40,jsamp=40 REAL(KIND=wp), DIMENSION(jpi,jpj) :: & & zerrrd ALLOCATE( & & ztem ( jpi, jpj ), & & zsal ( jpi, jpj ), & & zdep ( jpi, jpj ), & & zrd_out( jpi, jpj ), & & zrd_tl( jpi, jpj ), & & zs_tlin( jpi, jpj ), & & zt_tlin( jpi, jpj ), & & zrd_wop( jpi, jpj ), & & z2r ( jpi, jpj ) ) !-------------------------------------------------------------------- ! Reset the tangent and adjoint variables !-------------------------------------------------------------------- ztem ( :,:) = 0.0_wp zsal ( :,:) = 0.0_wp zt_tlin( :,:) = 0.0_wp zs_tlin( :,:) = 0.0_wp zrd_out( :,:) = 0.0_wp zrd_wop( :,:) = 0.0_wp zscerrrd( :) = 0.0_wp zscrd(:) = 0.0_wp IF ( tlm_bch == 2 ) zrd_tl ( :,:) = 0.0_wp !-------------------------------------------------------------------- ! Output filename Xn=F(X0) !-------------------------------------------------------------------- !! CALL tlm_namrd gamma = h_ratio file_wop='trj_wop_eos_2d' file_xdx='trj_xdx_eos_2d' !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN CALL grid_rd_sd( 596035, z2r, 'T', 0.0_wp, stdt) DO jj = nldj, nlej DO ji = nldi, nlei zt_tlin(ji,jj) = z2r(ji,jj) END DO END DO CALL grid_rd_sd( 371836, z2r, 'S', 0.0_wp, stds) DO jj = nldj, nlej DO ji = nldi, nlei zs_tlin(ji,jj) = z2r(ji,jj) END DO END DO ENDIF !-------------------------------------------------------------------- ! Complete Init for Direct !------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL istate_p ! *** initialize the reference trajectory ! ------------ CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) ! Initialize the reference state ztem(:,:) = tn(:,:,2) zsal(:,:) = sn(:,:,2) zdep(:,:) = fsdept(:,:,2) IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN zt_tlin(:,:) = gamma * zt_tlin(:,:) ztem(:,:) = ztem(:,:) + zt_tlin(:,:) zs_tlin(:,:) = gamma * zs_tlin(:,:) zsal(:,:) = zsal(:,:) + zs_tlin(:,:) ENDIF !-------------------------------------------------------------------- ! Compute the direct model F(X0,t=n) = Xn !-------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL eos(ztem, zsal, zdep, zrd_out) rhd (:,:,2) = zrd_out(:,:) 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: dy^* = W dy !-------------------------------------------------------------------- CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) ztem(:,:) = tn(:,:,2) zsal(:,:) = sn(:,:,2) !----------------------------------------------------------------------- ! Initialization of the dynamics and tracer fields for the tangent !----------------------------------------------------------------------- CALL eos_insitu_2d_tan(ztem, zsal, zdep, zt_tlin, zs_tlin, zrd_tl) !-------------------------------------------------------------------- ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) !-------------------------------------------------------------------- zsp2 = DOT_PRODUCT( zrd_tl, zrd_tl ) !-------------------------------------------------------------------- ! Storing data !-------------------------------------------------------------------- CALL trj_rd_spl(file_wop) zrd_wop (:,:) = rhd (:,:,2) CALL trj_rd_spl(file_xdx) zrd_out (:,:) = rhd (:,:,2) !-------------------------------------------------------------------- ! 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 jj = 1, jpj DO ji = 1, jpi zrd_out (ji,jj) = zrd_out (ji,jj) - zrd_wop (ji,jj) zrd_wop (ji,jj) = zrd_out (ji,jj) - zrd_tl (ji,jj) IF ( zrd_tl(ji,jj) .NE. 0.0_wp ) & & zerrrd(ji,jj) = zrd_out(ji,jj)/zrd_tl(ji,jj) IF( (MOD(ji, isamp) .EQ. 0) .AND. & & (MOD(jj, jsamp) .EQ. 0) ) THEN ii = ii+1 iiposrd(ii) = ji ijposrd(ii) = jj IF ( INT(tmask(ji,jj,2)) .NE. 0) THEN zscrd (ii) = zrd_wop(ji,jj) zscerrrd (ii) = ( zerrrd( ji,jj) - 1.0_wp ) / gamma ENDIF ENDIF END DO END DO zsp1 = DOT_PRODUCT( zrd_out, zrd_out ) zsp3 = DOT_PRODUCT( zrd_wop, zrd_wop ) !-------------------------------------------------------------------- ! Print the linearization error En - norme 2 !-------------------------------------------------------------------- ! 14 char:'12345678901234' cl_name = 'eos_2d:En ' zzsp = SQRT(zsp3) zgsp5 = zzsp CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Compute TLM norm2 !-------------------------------------------------------------------- zzsp = SQRT(zsp2) zgsp4 = zzsp cl_name = 'eos_2d:Ln2' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Print the linearization error Nn - norme 2 !-------------------------------------------------------------------- zzsp = SQRT(zsp1) cl_name = 'eos_2d: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) 'eosins2d', 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 = 'eosins2d' IF(lwp) THEN DO ii=1, 100, 1 IF ( zscrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrd ', & & cur_loop, h_ratio, ii, iiposrd(ii), ijposrd(ii), zscrd(ii) ENDDO DO ii=1, 100, 1 IF ( zscerrrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrrd', & & cur_loop, h_ratio, ii, iiposrd(ii), ijposrd(ii), zscerrrd(ii) ENDDO ! write separator WRITE(numtan_sc,"(A4)") '====' ENDIF ENDIF DEALLOCATE( & & zrd_out, zrd_tl, zrd_wop, & & ztem, zsal, zdep, & & zt_tlin, zs_tlin, z2r & & ) END SUBROUTINE eos_insitu_2d_tlm_tst SUBROUTINE eos_pot_1pt_tlm_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_pot_1pt_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 eosbn2, ONLY: & ! horizontal & vertical advective trend & eos USE tamtrj ! writing out state trajectory USE par_tlm, ONLY: & & tlm_bch, & & cur_loop, & & h_ratio USE istate_mod USE trj_tam USE oce , ONLY: & ! ocean dynamics and tracers variables & tn, sn, rhd, rhop USE in_out_manager, ONLY: & ! I/O manager & nitend, & & nit000 USE tamctl, ONLY: & ! Control parameters & numtan, numtan_sc !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit INTEGER:: & & numsctlm !! * Local declarations REAL(KIND=wp) :: & & ztem, & & zsal, & & zt_tlin, & & zs_tlin, & & zrh_out, & & zrh_tl, & & zrh_wop REAL(KIND=wp) :: & & zsp1, zsp2, & & zsp3, zzsp, & & gamma, & & zgsp1, & & zgsp2, & & zgsp3, & & zgsp4, & & zgsp5, & & zgsp6, & & zgsp7 CHARACTER(LEN=14) :: cl_name CHARACTER (LEN=128) :: file_out, file_wop, file_xdx CHARACTER (LEN=90) :: FMT ! Initialize the reference state ztem = 23.7 zsal = 30.1 !-------------------------------------------------------------------- ! Reset the tangentvariables !-------------------------------------------------------------------- zt_tlin = 1.12_wp zs_tlin = 0.123_wp zrh_out = 0.0_wp zrh_wop = 0.0_wp IF ( tlm_bch == 2 ) zrh_tl = 0.0_wp !-------------------------------------------------------------------- ! Output filename Xn=F(X0) !-------------------------------------------------------------------- !! CALL tlm_namrd gamma = h_ratio file_wop='trj_wop_eos_1pt' file_xdx='trj_xdx_eos_1pt' !-------------------------------------------------------------------- ! Complete Init for Direct !------------------------------------------------------------------- IF ( tlm_bch /= 2 ) 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 zt_tlin = gamma * zt_tlin ztem = ztem + zt_tlin zs_tlin = gamma * zs_tlin zsal = zsal + zs_tlin ENDIF !-------------------------------------------------------------------- ! Compute the direct model F(X0,t=n) = Xn !-------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL eos(ztem, zsal, zrh_out) rhop (1,1,1) = zrh_out 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: dy^* = W dy !-------------------------------------------------------------------- CALL trj_rea( nit000-1, 1 ) CALL trj_rea( nit000, 1 ) ztem = 23.7 zsal = 30.1 !----------------------------------------------------------------------- ! Initialization of the dynamics and tracer fields for the tangent !----------------------------------------------------------------------- CALL eos_pot_1pt_tan(ztem, zsal, zt_tlin, zs_tlin, zrh_tl) !-------------------------------------------------------------------- ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) !-------------------------------------------------------------------- zsp2 = abs(zrh_tl) !-------------------------------------------------------------------- ! Storing data !-------------------------------------------------------------------- CALL trj_rd_spl(file_wop) zrh_wop = rhop (1,1,1) CALL trj_rd_spl(file_xdx) zrh_out = rhop (1,1,1) !-------------------------------------------------------------------- ! 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 zrh_out = zrh_out - zrh_wop zrh_wop = zrh_out - zrh_tl zsp1 = abs(zrh_out) zsp3 = abs(zrh_wop) !-------------------------------------------------------------------- ! Print the linearization error En - norme 2 !-------------------------------------------------------------------- ! 14 char:'12345678901234' cl_name = 'eos_1pt:En ' zzsp = zsp3 zgsp5 = zzsp CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Compute TLM norm2 !-------------------------------------------------------------------- zzsp = zsp2 zgsp4 = zzsp cl_name = 'eos_1pt:Ln2' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Print the linearization error Nn - norme 2 !-------------------------------------------------------------------- zzsp = zsp1 cl_name = 'eos_1pt: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) 'eos1pt ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7 ENDIF END SUBROUTINE eos_pot_1pt_tlm_tst SUBROUTINE bn2_tlm_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE bn2_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-12 (A. Vigilant) !!----------------------------------------------------------------------- !! * Modules used USE eosbn2, ONLY: & ! horizontal & vertical advective trend & bn2 USE tamtrj ! writing out state trajectory USE par_tlm, ONLY: & & tlm_bch, & & cur_loop, & & h_ratio USE istate_mod USE sshwzv ! vertical velocity USE gridrandom, ONLY: & & grid_rd_sd USE trj_tam USE oce , ONLY: & ! ocean dynamics and tracers variables & tn, sn, rn2 USE oce_tam , ONLY: & & tn_tl, & & sn_tl, & & rn2_tl USE in_out_manager, ONLY: & ! I/O manager & & nit000 USE tamctl, ONLY: & ! Control parameters & numtan, numtan_sc !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit !! * Local declarations INTEGER :: & & ji, & ! dummy loop indices & jj, & & jk REAL(KIND=wp) :: & & zsp1, & ! scalar product involving the tangent routine & zsp2, & ! scalar product involving the tangent routine & zsp3, & ! scalar product involving the tangent routine & zzsp, & ! scalar product involving the tangent routine & gamma, & & zgsp1, & & zgsp2, & & zgsp3, & & zgsp4, & & zgsp5, & & zgsp6, & & zgsp7 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & & ztn_tlin, & ! potential temperature & zsn_tlin, & ! salinity & zrn2_out, & ! Brunt-Vaisala frequency [s-1] Direct output & zrn2_wop, & ! Brunt-Vaisala frequency [s-1] Direct output w/o perturbation & z3r CHARACTER(LEN=14) :: cl_name CHARACTER (LEN=128) :: file_out, file_wop, file_xdx CHARACTER (LEN=90) :: FMT REAL(KIND=wp), DIMENSION(100):: & & zscrn2, & & zscerrrn2 INTEGER, DIMENSION(100):: & & iiposrn2,ijposrn2, ikposrn2 INTEGER:: & & ii, & & isamp=40, & & jsamp=40, & & ksamp=10, & & numsctlm REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: & & zerrrn2 ! Allocate memory ALLOCATE( & & zsn_tlin( jpi, jpj, jpk ), & & ztn_tlin( jpi, jpj, jpk ), & & zrn2_out( jpi, jpj, jpk ), & & zrn2_wop( jpi, jpj, jpk ), & & z3r( jpi, jpj, jpk ) ) !-------------------------------------------------------------------- ! Output filename Xn=F(X0) !-------------------------------------------------------------------- !! CALL tlm_namrd gamma = h_ratio file_wop='trj_wop_bn2' file_xdx='trj_xdx_bn2' !-------------------------------------------------------------------- ! Initialize the tangent input with random noise: dx !-------------------------------------------------------------------- IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN CALL grid_rd_sd( 264940, z3r, 'T', 0.0_wp, stdt) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO CALL grid_rd_sd( 618304, z3r, 'T', 0.0_wp, stds) DO jk = 1, jpk DO jj = nldj, nlej DO ji = nldi, nlei zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk) END DO END DO END DO ENDIF !-------------------------------------------------------------------- ! Complete Init for Direct !------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL istate_p ! *** initialize the reference trajectory ! ------------ CALL trj_rea( nit000-1, 1 ) IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN ztn_tlin(:,:,:) = gamma * ztn_tlin(:,:,:) tn(:,:,:) = tn(:,:,:) + ztn_tlin(:,:,:) zsn_tlin(:,:,:) = gamma * zsn_tlin(:,:,:) sn(:,:,:) = sn(:,:,:) + zsn_tlin(:,:,:) ENDIF !-------------------------------------------------------------------- ! Compute the direct model F(X0,t=n) = Xn !-------------------------------------------------------------------- IF ( tlm_bch /= 2 ) CALL bn2(tn, sn, rn2) 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: dy^* = W dy !-------------------------------------------------------------------- CALL trj_rea( nit000-1, 1 ) tn_tl (:,:,:) = ztn_tlin (:,:,:) sn_tl (:,:,:) = zsn_tlin (:,:,:) !----------------------------------------------------------------------- ! Initialization of the dynamics and tracer fields for the tangent !----------------------------------------------------------------------- CALL bn2_tan(tn, sn, ztn_tlin, zsn_tlin, rn2_tl) !-------------------------------------------------------------------- ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) !-------------------------------------------------------------------- zsp2 = DOT_PRODUCT( rn2_tl, rn2_tl ) !-------------------------------------------------------------------- ! Storing data !-------------------------------------------------------------------- CALL trj_rd_spl(file_wop) zrn2_wop (:,:,:) = rn2 (:,:,:) CALL trj_rd_spl(file_xdx) zrn2_out (:,:,:) = rn2 (:,:,:) !-------------------------------------------------------------------- ! 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 zrn2_out (ji,jj,jk) = zrn2_out (ji,jj,jk) - zrn2_wop (ji,jj,jk) zrn2_wop (ji,jj,jk) = zrn2_out (ji,jj,jk) - rn2_tl (ji,jj,jk) IF ( rn2_tl(ji,jj,jk) .NE. 0.0_wp ) & & zerrrn2(ji,jj,jk) = zrn2_out(ji,jj,jk)/rn2_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 iiposrn2(ii) = ji ijposrn2(ii) = jj ikposrn2(ii) = jk IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN zscrn2 (ii) = zrn2_wop(ji,jj,jk) zscerrrn2 (ii) = ( zerrrn2(ji,jj,jk) - 1.0_wp ) / gamma ENDIF ENDIF END DO END DO END DO zsp1 = DOT_PRODUCT( zrn2_out, zrn2_out ) zsp3 = DOT_PRODUCT( zrn2_wop, zrn2_wop ) !-------------------------------------------------------------------- ! Print the linearization error En - norme 2 !-------------------------------------------------------------------- ! 14 char:'12345678901234' cl_name = 'eosbn2_tam:En ' zzsp = SQRT(zsp3) zgsp5 = zzsp CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Compute TLM norm2 !-------------------------------------------------------------------- zzsp = SQRT(zsp2) zgsp4 = zzsp cl_name = 'eosbn2_tam:Ln2' CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) !-------------------------------------------------------------------- ! Print the linearization error Nn - norme 2 !-------------------------------------------------------------------- zzsp = SQRT(zsp1) cl_name = 'traadv: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) 'eosbn2 ', 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 = 'eosbn2 ' IF(lwp) THEN DO ii=1, 100, 1 IF ( zscrn2(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrn2 ', & & cur_loop, h_ratio, ii, iiposrn2(ii), ijposrn2(ii), zscrn2(ii) ENDDO ! write separator WRITE(numtan_sc,"(A4)") '====' ENDIF ENDIF DEALLOCATE( & & ztn_tlin, zsn_tlin, & & zrn2_out, zrn2_wop, & & z3r & & ) END SUBROUTINE bn2_tlm_tst SUBROUTINE eos_tlm_tst( kumadt ) !!----------------------------------------------------------------------- !! !! *** ROUTINE eos_tlm_tst *** !! !! ** Purpose : Test the tangent routine. !! !! History : !! ! 09-08 (F. Vigilant) !!----------------------------------------------------------------------- !! * Modules used !! * Arguments INTEGER, INTENT(IN) :: & & kumadt ! Output unit CALL eos_insitu_tlm_tst( kumadt ) CALL eos_insitu_pot_tlm_tst( kumadt ) CALL eos_insitu_2d_tlm_tst( kumadt ) CALL eos_pot_1pt_tlm_tst( kumadt ) END SUBROUTINE eos_tlm_tst #endif #endif !!====================================================================== END MODULE eosbn2_tam