New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
eosbn2_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA/eosbn2_tam.F90 @ 3640

Last change on this file since 3640 was 3640, checked in by pabouttier, 10 years ago

Missing allocation/deallocation in TAM routines - See ticket #1013

  • Property svn:executable set to *
File size: 140.3 KB
Line 
1MODULE eosbn2_tam
2   !!==============================================================================
3   !!                       ***  MODULE  eosbn2_tam  ***
4   !! Ocean diagnostic variable : equation of state - in situ and potential density
5   !!                                               - Brunt-Vaisala frequency
6   !!                            Tangent and Adjoint module
7   !!==============================================================================
8   !! History :  OPA  ! 1989-03  (O. Marti)  Original code
9   !!            6.0  ! 1994-07  (G. Madec, M. Imbard)  add bn2
10   !!            6.0  ! 1994-08  (G. Madec)  Add Jackett & McDougall eos
11   !!            7.0  ! 1996-01  (G. Madec)  statement function for e3
12   !!            8.1  ! 1997-07  (G. Madec)  density instead of volumic mass
13   !!             -   ! 1999-02  (G. Madec, N. Grima) semi-implicit pressure gradient
14   !!            8.2  ! 2001-09  (M. Ben Jelloul)  bugfix on linear eos
15   !!   NEMO     1.0  ! 2002-10  (G. Madec)  add eos_init
16   !!             -   ! 2002-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d
17   !!             -   ! 2003-08  (G. Madec)  F90, free form
18   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function
19   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
20   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp
21   !!   History of the TAM Module :
22   !!             8.2 ! 2005-03  (F. Van den Berghe, A. Weaver, N. Daget)
23   !!                                              - eostan.F
24   !!             9.0 ! 2007-07  (K. Mogensen) Initial version based on eostan.F
25   !!                 ! 2008-07 (A. Vidard) bug fix in computation of prd_tl if neos=1
26   !!                 ! 2008-11  (A. Vidard) TAM of the 06-08 version
27   !!  NEMO       3.2 ! 2010-04  (F. Vigilant) version 3.2
28   !!             3.4 ! 2012-04  (P.-A. Bouttier) version 3.4
29   !!----------------------------------------------------------------------
30
31   !!----------------------------------------------------------------------
32   !!  Direct subroutines
33   !!   eos            : generic interface of the equation of state
34   !!   eos_insitu     : Compute the in situ density
35   !!   eos_insitu_pot : Compute the insitu and surface referenced potential
36   !!                    volumic mass
37   !!   eos_insitu_2d  : Compute the in situ density for 2d fields
38   !!   eos_bn2        : Compute the Brunt-Vaisala frequency
39   !!   eos_alpbet     : calculates the in situ thermal/haline expansion ratio
40   !!   tfreez         : Compute the surface freezing temperature
41   !!   eos_init       : set eos parameters (namelist)
42   !!----------------------------------------------------------------------
43   !! * Modules used
44#if defined key_zdfddm
45   USE oce_tam
46#endif
47   USE dom_oce
48   USE par_kind
49   USE par_oce
50   USE oce
51   USE phycst
52   USE in_out_manager
53#if defined key_zdfddm
54   USE zdfddm          ! vertical physics: double diffusion
55#endif
56   USE eosbn2
57   USE gridrandom
58   USE dotprodfld
59   USE tstool_tam
60   USE wrk_nemo
61   USE timing
62   USE lib_mpp
63   IMPLICIT NONE
64   PRIVATE
65   !
66   !! * Interface
67    INTERFACE eos_tan
68       MODULE PROCEDURE eos_insitu_tan, eos_insitu_pot_tan, eos_insitu_2d_tan, &
69                        eos_alpbet_tan
70    END INTERFACE
71    INTERFACE eos_adj
72       MODULE PROCEDURE eos_insitu_adj, eos_insitu_pot_adj, eos_insitu_2d_adj, &
73                        eos_alpbet_adj
74    END INTERFACE
75   INTERFACE bn2_tan
76      MODULE PROCEDURE eos_bn2_tan
77   END INTERFACE
78   INTERFACE bn2_adj
79      MODULE PROCEDURE eos_bn2_adj
80   END INTERFACE
81   !
82   !! * Routine accessibility
83   PUBLIC eos_tan    ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90
84   PUBLIC bn2_tan    ! called by step.F90
85   PUBLIC eos_adj    ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90
86   PUBLIC bn2_adj    ! called by step.F90
87#if defined key_tam
88   PUBLIC eos_adj_tst
89   PUBLIC bn2_adj_tst
90#endif
91   !! * Substitutions
92#  include "domzgr_substitute.h90"
93#  include "vectopt_loop_substitute.h90"
94   !!----------------------------------------------------------------------
95   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
96   !! $Id$
97   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
98   !!----------------------------------------------------------------------
99CONTAINS
100   SUBROUTINE eos_insitu_tan( pts, pts_tl, prd_tl )
101      !!-----------------------------------------------------------------------
102      !!
103      !!        ***  ROUTINE eos_insitu_tan : TL OF ROUTINE   eos_insitu ***
104      !!
105      !! ** Purpose of direct routine   : Compute the in situ density
106      !!       (ratio rho/rau0) from potential temperature and salinity
107      !!       using an equation of state defined through the namelist
108      !!       parameter nn_eos.
109      !!
110      !! ** Method of direct routine    : 3 cases:
111      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
112      !!         the in situ density is computed directly as a function of
113      !!         potential temperature relative to the surface (the opa t
114      !!         variable), salt and pressure (assuming no pressure variation
115      !!         along geopotential surfaces, i.e. the pressure p in decibars
116      !!         is approximated by the depth in meters.
117      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
118      !!         with pressure                      p        decibars
119      !!              potential temperature         t        deg celsius
120      !!              salinity                      s        psu
121      !!              reference volumic mass        rau0     kg/m**3
122      !!              in situ volumic mass          rho      kg/m**3
123      !!              in situ density anomalie      prd      no units
124      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
125      !!          t = 40 deg celcius, s=40 psu
126      !!      nn_eos = 1 : linear equation of state function of temperature only
127      !!              prd(t) = 0.0285 - rn_alpha * t
128      !!      nn_eos = 2 : linear equation of state function of temperature and
129      !!               salinity
130      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
131      !!      Note that no boundary condition problem occurs in this routine
132      !!      as (ptem,psal) are defined over the whole domain.
133      !!
134      !! ** Comments on Adjoint Routine :
135      !!      Care has been taken to avoid division by zero when computing
136      !!      the inverse of the square root of salinity at masked salinity
137      !!      points.
138      !!
139      !! * Arguments
140      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts, &   ! 1 : potential temperature  [Celcius]
141      !                                                         ! 2 : salinity               [psu]
142      &                                                pts_tl   ! 1 : TL of potential temperature [Celsius]
143                                                                ! 2 : TL of salinity [psu]
144      REAL(wp), DIMENSION(:,:,:), INTENT( out ) ::   &
145         & prd_tl                   ! TL of potential density (surface referenced)
146      !! * Local declarations
147      INTEGER ::  ji, jj, jk      ! dummy loop indices
148      REAL(wp) ::   &             ! temporary scalars
149         zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw,   &
150         zb, zd, zc, zaw, za, zb1, za1, zkw, zk0,               &
151         zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl,          &
152         zr4tl, zrhoptl, zetl, zbwtl,                           &
153         zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl,           &
154         zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps,                &
155         zmask, zrau0r
156      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
157      !!----------------------------------------------------------------------
158      !
159      IF( nn_timing == 1 ) CALL timing_start('eos_tan')
160      !
161      CALL wrk_alloc( jpi, jpj, jpk, zws )
162      !
163      SELECT CASE ( nn_eos )
164
165      CASE ( 0 )               !== Jackett and McDougall (1994) formulation ==!
166         zrau0r = 1._wp / rau0
167#ifdef key_sp
168         zeps = 1.e-7_wp
169#else
170         zeps = 1.e-14_wp
171#endif
172!CDIR NOVERRCHK
173         zws(:,:,:) = SQRT( ABS( pts(:,:,:, jp_sal) ) )
174         !
175         DO jk = 1, jpkm1
176            DO jj = 1, jpj
177               DO ji = 1, jpi
178                  zt = pts(ji,jj,jk,jp_tem)
179                  zs = pts(ji,jj,jk,jp_sal)
180                  zh = fsdept(ji,jj,jk)        ! depth
181                  zsr= zws(ji,jj,jk)           ! square root salinity
182                  ! compute volumic mass pure water at atm pressure
183                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp)*zt   &
184                     -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
185                  ! seawater volumic mass atm pressure
186                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   &
187                     -4.0899e-3_wp ) *zt+0.824493_wp
188                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp
189                  zr4= 4.8314e-4_wp
190                  !
191                  ! potential volumic mass (reference to the surface)
192                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
193                  !
194                  ! add the compression terms
195                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
196                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
197                  zb = zbw + ze * zs
198
199                  zd = -2.042967e-2_wp
200                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
201                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
202                  za = ( zd*zsr + zc ) *zs + zaw
203
204                  zb1=   (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp
205                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)*zt-65.00517_wp ) *zt+1044.077_wp
206                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
207                  zk0= ( zb1*zsr + za1 )*zs + zkw
208
209                  ! Tangent linear part
210
211                  zttl = pts_tl(ji,jj,jk, jp_tem)
212                  zstl = pts_tl(ji,jj,jk, jp_sal)
213
214                  zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
215                     &   * tmask(ji,jj,jk)                * zstl
216
217                  zr1tl= ( ( ( (  5.*6.536332e-9_wp   * zt &
218                     &           -4.*1.120083e-6_wp ) * zt &
219                     &           +3.*1.001685e-4_wp ) * zt &
220                     &           -2.*9.095290e-3_wp ) * zt &
221                     &           +   6.793952e-2_wp ) * zttl
222
223                  zr2tl= ( ( (    4.*5.3875e-9_wp   * zt &
224                     &           -3.*8.2467e-7_wp ) * zt &
225                     &           +2.*7.6438e-5_wp ) * zt &
226                     &           -   4.0899e-3_wp ) * zttl
227
228                  zr3tl= (       -2.*1.6546e-6_wp   * zt &
229                     &           +   1.0227e-4_wp ) * zttl
230
231                  zrhoptl=                                  zr1tl &
232                     &     + zs                           * zr2tl &
233                     &     + zsr * zs                     * zr3tl &
234                     &     + zr3 * zs                     * zsrtl &
235                     &     + (  2. * zr4 * zs + zr2               &
236                     &        + zr3 * zsr                         )     * zstl
237
238                  zetl = (       -2.*3.508914e-8_wp   * zt &
239                     &           -   1.248266e-8_wp        ) * zttl
240
241                  zbwtl= (        2.*1.296821e-6_wp   * zt &
242                     &           -   5.782165e-9_wp        ) * zttl
243
244                  zbtl =                                    zbwtl &
245                     &    + zs                            * zetl  &
246                  &       + ze                            * zstl
247
248                  zctl = (       -2.*7.267926e-5_wp   * zt &
249                     &           +   2.598241e-3_wp        ) * zttl
250
251                  zawtl= ( (      3.*5.939910e-6_wp   * zt &
252                     &           +2.*2.512549e-3_wp ) * zt &
253                     &           -   0.1028859_wp          ) * zttl
254
255                  zatl =                                    zawtl &
256                  &      + zd * zs                        * zsrtl &
257                  &      + zs                             * zctl  &
258                  &      + ( zd * zsr + zc )              * zstl
259
260                  zb1tl= (       -2.*0.1909078_wp     * zt &
261                     &           +   7.390729_wp           ) * zttl
262
263                  za1tl= ( (      3.*2.326469e-3_wp   * zt &
264                     &           +2.*1.553190_wp    ) * zt &
265                     &           -   65.00517_wp           ) * zttl
266
267                  zkwtl= ( ( (   -4.*1.361629e-4_wp   * zt &
268                     &           -3.*1.852732e-2_wp ) * zt &
269                     &           -2.*30.41638_wp    ) * zt &
270                     &           +   2098.925_wp           ) * zttl
271
272                  zk0tl=                                    zkwtl &
273                     &  + zb1 * zs                        * zsrtl &
274                     &  + zs  * zsr                       * zb1tl &
275                     &  + zs                              * za1tl &
276                     &  + ( zb1 * zsr + za1 )             * zstl
277
278                  ! Masked in situ density anomaly
279
280                  zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
281                  zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
282
283                  prd_tl(ji,jj,jk) = tmask(ji,jj,jk) * zrdc2 * &
284                     &               (                              zrhoptl &
285                     &                 - zrdc2 * zh * zrdc1**2 * zrhop      &
286                     &                   * (                        zk0tl   &
287                     &                       - zh * (               zatl    &
288                     &                                - zh        * zbtl ) ) )&
289                     &               * zrau0r
290               END DO
291            END DO
292         END DO
293         !
294      CASE ( 1 )               !== Linear formulation function of temperature only ==!
295         DO jk = 1, jpkm1
296            prd_tl(:,:,jk) = ( - rn_alpha * pts_tl(:,:,jk,jp_tem) ) * tmask(:,:,jk)
297         END DO
298         !
299      CASE ( 2 )               !== Linear formulation function of temperature and salinity ==!
300
301         !                                                ! ===============
302         DO jk = 1, jpkm1
303            prd_tl(:,:,jk) = ( rn_beta  * pts_tl(:,:,jk,jp_sal) - rn_alpha * pts_tl(:,:,jk,jp_tem ) ) * tmask(:,:,jk)
304         END DO
305         !
306      END SELECT
307      !
308      CALL wrk_dealloc( jpi, jpj, jpk, zws )
309      !
310      IF( nn_timing == 1 ) CALL timing_stop('eos_tan')
311      !
312   END SUBROUTINE eos_insitu_tan
313
314   SUBROUTINE eos_insitu_pot_tan( pts, pts_tl, prd_tl, prhop_tl)
315      !!----------------------------------------------------------------------
316      !!                  ***  ROUTINE eos_insitu_pot_tan  ***
317      !!
318      !! ** Purpose or the direct routine:
319      !!      Compute the in situ density (ratio rho/rau0) and the
320      !!      potential volumic mass (Kg/m3) from potential temperature and
321      !!      salinity fields using an equation of state defined through the
322      !!     namelist parameter nn_eos.
323      !!
324      !! ** Method  :
325      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
326      !!         the in situ density is computed directly as a function of
327      !!         potential temperature relative to the surface (the opa t
328      !!         variable), salt and pressure (assuming no pressure variation
329      !!         along geopotential surfaces, i.e. the pressure p in decibars
330      !!         is approximated by the depth in meters.
331      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
332      !!              rhop(t,s)  = rho(t,s,0)
333      !!         with pressure                      p        decibars
334      !!              potential temperature         t        deg celsius
335      !!              salinity                      s        psu
336      !!              reference volumic mass        rau0     kg/m**3
337      !!              in situ volumic mass          rho      kg/m**3
338      !!              in situ density anomalie      prd      no units
339      !!
340      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
341      !!          t = 40 deg celcius, s=40 psu
342      !!
343      !!      nn_eos = 1 : linear equation of state function of temperature only
344      !!              prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t
345      !!              rhop(t,s)  = rho(t,s)
346      !!
347      !!      nn_eos = 2 : linear equation of state function of temperature and
348      !!               salinity
349      !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0
350      !!                       = rn_beta * s - rn_alpha * tn - 1.
351      !!              rhop(t,s)  = rho(t,s)
352      !!      Note that no boundary condition problem occurs in this routine
353      !!      as (tn,sn) or (ta,sa) are defined over the whole domain.
354      !!
355      !! ** Action  : - prd  , the in situ density (no units)
356      !!              - prhop, the potential volumic mass (Kg/m3)
357      !!
358      !! References :
359      !!      Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
360      !!      Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
361      !!
362      !!----------------------------------------------------------------------
363      !! * Arguments
364      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts, &  ! 1 : potential temperature  [Celcius]
365      !                                                                 ! 2 : salinity               [psu]
366      &                                                         pts_tl  ! 1 : TL of potential temperature  [Celcius]
367      !                                                                 ! 2 : TL of salinity               [psu]
368      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd_tl  ! TL of in_situ density [-]
369      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop_tl  ! TL of potential density (surface referenced)
370      !! * Local declarations
371      INTEGER ::  ji, jj, jk      ! dummy loop indices
372      REAL(wp) ::   &             ! temporary scalars
373         zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw,   &
374         zb, zd, zc, zaw, za, zb1, za1, zkw, zk0,               &
375         zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl,          &
376         zr4tl, zrhoptl, zetl, zbwtl,                           &
377         zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl,           &
378         zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps,                &
379         zmask, zrau0r
380      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
381      !!----------------------------------------------------------------------
382      !
383      IF( nn_timing == 1 ) CALL timing_start('eos-p_tan')
384      !
385      CALL wrk_alloc( jpi, jpj, jpk, zws )
386      !
387      SELECT CASE ( nn_eos )
388      !
389      CASE ( 0 )               !== Jackett and McDougall (1994) formulation ==!
390         zrau0r = 1.e0 / rau0
391#ifdef key_sp
392         zeps = 1.e-7
393#else
394         zeps = 1.e-14
395#endif
396!CDIR NOVERRCHK
397         zws(:,:,:) = SQRT( ABS( pts(:,:,:, jp_sal) ) )
398         !
399         DO jk = 1, jpkm1
400            DO jj = 1, jpj
401               DO ji = 1, jpi
402                  zt = pts(ji,jj,jk, jp_tem)
403                  zs = pts(ji,jj,jk, jp_sal)
404                  zh = fsdept(ji,jj,jk)        ! depth
405                  zsr = zws(ji,jj,jk)          ! square root salinity
406                  ! compute volumic mass pure water at atm pressure
407                  zr1 = ( ( ( ( 6.536332e-9_wp   * zt - 1.120083e-6_wp ) * zt   &
408                     &        + 1.001685e-4_wp ) * zt - 9.095290e-3_wp ) * zt   &
409                     &        + 6.793952e-2_wp ) * zt + 999.842594_wp
410                  ! seawater volumic mass atm pressure
411                  zr2 = ( ( ( 5.3875e-9_wp   * zt - 8.2467e-7_wp ) * zt         &
412                     &      + 7.6438e-5_wp ) * zt - 4.0899e-3_wp ) * zt         &
413                     &      + 0.824493_wp
414                  zr3 = ( -1.6546e-6_wp * zt + 1.0227e-4_wp ) * zt - 5.72466e-3_wp
415                  zr4 = 4.8314e-4_wp
416
417                  ! potential volumic mass (reference to the surface)
418                  zrhop= ( zr4 * zs + zr3 * zsr + zr2 ) * zs + zr1
419
420                  ! add the compression terms
421                  ze  = ( -3.508914e-8_wp * zt - 1.248266e-8_wp ) * zt - 2.595994e-6_wp
422                  zbw = (  1.296821e-6_wp * zt - 5.782165e-9_wp ) * zt + 1.045941e-4_wp
423                  zb  = zbw + ze * zs
424
425                  zd = -2.042967e-2_wp
426                  zc =   (-7.267926e-5_wp * zt + 2.598241e-3_wp ) * zt + 0.1571896_wp
427                  zaw= ( ( 5.939910e-6_wp * zt + 2.512549e-3_wp ) * zt - 0.1028859_wp ) * zt - 4.721788_wp
428                  za = ( zd * zsr + zc ) * zs + zaw
429
430                  zb1 =   (-0.1909078_wp   * zt + 7.390729_wp ) * zt - 55.87545_wp
431                  za1 = ( ( 2.326469e-3_wp * zt + 1.553190_wp ) * zt - 65.00517_wp  &
432                     &  ) * zt + 1044.077_wp
433                  zkw = ( ( (-1.361629e-4_wp * zt - 1.852732e-2_wp ) * zt - 30.41638_wp &
434                     &    ) * zt + 2098.925_wp ) * zt + 190925.6_wp
435                  zk0 = ( zb1 * zsr + za1 ) * zs + zkw
436
437
438                  zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
439                  zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
440
441                  ! Tangent linear part
442
443                  zttl = pts_tl(ji,jj,jk, jp_tem)
444                  zstl = pts_tl(ji,jj,jk, jp_sal)
445
446                  zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
447                     &   * tmask(ji,jj,jk)                * zstl
448
449                  zr1tl= ( ( ( (  5.*6.536332e-9_wp   * zt &
450                     &           -4.*1.120083e-6_wp ) * zt &
451                     &           +3.*1.001685e-4_wp ) * zt &
452                     &           -2.*9.095290e-3_wp ) * zt &
453                     &           +   6.793952e-2_wp        ) * zttl
454
455                  zr2tl= ( ( (    4.*5.3875e-9_wp   * zt &
456                     &           -3.*8.2467e-7_wp ) * zt &
457                     &           +2.*7.6438e-5_wp ) * zt &
458                     &           -   4.0899e-3_wp        ) * zttl
459
460                  zr3tl= (       -2.*1.6546e-6_wp   * zt &
461                     &           +   1.0227e-4_wp        ) * zttl
462
463                  zrhoptl=                                  zr1tl &
464                     &     + zs                           * zr2tl &
465                     &     + zsr * zs                     * zr3tl &
466                     &     + zr3 * zs                     * zsrtl &
467                     &     + (  2. * zr4 * zs + zr2               &
468                     &        + zr3 * zsr           )     * zstl
469
470                  prhop_tl(ji,jj,jk) = zrhoptl * tmask(ji,jj,jk)
471
472                  zetl = (       -2.*3.508914e-8_wp   * zt &
473                     &           -   1.248266e-8_wp        ) * zttl
474
475                  zbwtl= (        2.*1.296821e-6_wp   * zt &
476                     &           -   5.782165e-9_wp        ) * zttl
477
478                  zbtl =                                    zbwtl &
479                     &    + zs                            * zetl  &
480                  &       + ze                            * zstl
481
482                  zctl = (       -2.*7.267926e-5_wp   * zt &
483                     &           +   2.598241e-3_wp        ) * zttl
484
485                  zawtl= ( (      3.*5.939910e-6_wp   * zt &
486                     &           +2.*2.512549e-3_wp ) * zt &
487                     &           -   0.1028859_wp          ) * zttl
488
489                  zatl =                                    zawtl &
490                  &      + zd * zs                        * zsrtl &
491                  &      + zs                             * zctl  &
492                  &      + ( zd * zsr + zc )              * zstl
493
494                  zb1tl= (       -2.*0.1909078_wp     * zt &
495                     &           +   7.390729_wp           ) * zttl
496
497                  za1tl= ( (      3.*2.326469e-3_wp   * zt &
498                     &           +2.*1.553190_wp    ) * zt &
499                     &           -   65.00517_wp           ) * zttl
500
501                  zkwtl= ( ( (   -4.*1.361629e-4_wp   * zt &
502                     &           -3.*1.852732e-2_wp ) * zt &
503                     &           -2.*30.41638_wp    ) * zt &
504                     &           +   2098.925_wp           ) * zttl
505
506                  zk0tl=                                    zkwtl &
507                     &  + zb1 * zs                        * zsrtl &
508                     &  + zs  * zsr                       * zb1tl &
509                     &  + zs                              * za1tl &
510                     &  + ( zb1 * zsr + za1 )             * zstl
511
512                  ! Masked in situ density anomaly
513
514                  prd_tl(ji,jj,jk) = tmask(ji,jj,jk) * zrdc2 * &
515                     &               (                              zrhoptl &
516                     &                 - zrdc2 * zh * zrdc1**2 * zrhop      &
517                     &                   * (                        zk0tl   &
518                     &                       - zh * (               zatl    &
519                     &                                - zh        * zbtl ) ) )&
520                     &               * zrau0r
521               END DO
522            END DO
523         END DO
524         !
525      CASE ( 1 )                !==  Linear formulation = F( temperature )  ==!
526         DO jk = 1, jpkm1
527            prd_tl  (:,:,jk) = ( - rn_alpha * pts_tl(:,:,jk, jp_tem) ) * tmask(:,:,jk)
528            prhop_tl(:,:,jk) = ( rau0 * prd_tl(:,:,jk)        ) * tmask(:,:,jk)
529         END DO
530         !
531      CASE ( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
532         DO jk = 1, jpkm1
533            prd_tl(:,:,jk)   = ( rn_beta  * pts_tl(:,:,jk, jp_sal) - rn_alpha * pts_tl(:,:,jk, jp_tem) ) * tmask(:,:,jk)
534            prhop_tl(:,:,jk) = ( rau0 * prd_tl(:,:,jk) )       * tmask(:,:,jk)
535         END DO
536         !
537      END SELECT
538      !
539      CALL wrk_dealloc( jpi, jpj, jpk, zws )
540      !
541      IF( nn_timing == 1 ) CALL timing_stop('eos-p_tan')
542      !
543   END SUBROUTINE eos_insitu_pot_tan
544   SUBROUTINE eos_insitu_2d_tan( pts, pdep, pts_tl, prd_tl )
545      !!-----------------------------------------------------------------------
546      !!
547      !!        ***  ROUTINE eos_insitu_2d_tan : TL OF ROUTINE eos_insitu_2d ***
548      !!
549      !! ** Purpose of direct routine   : Compute the in situ density
550      !!       (ratio rho/rau0) from potential temperature and salinity
551      !!       using an equation of state defined through the namelist
552      !!       parameter nn_eos. * 2D field case
553      !!
554      !! ** Method of direct routine    : 3 cases:
555      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
556      !!         the in situ density is computed directly as a function of
557      !!         potential temperature relative to the surface (the opa t
558      !!         variable), salt and pressure (assuming no pressure variation
559      !!         along geopotential surfaces, i.e. the pressure p in decibars
560      !!         is approximated by the depth in meters.
561      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
562      !!         with pressure                      p        decibars
563      !!              potential temperature         t        deg celsius
564      !!              salinity                      s        psu
565      !!              reference volumic mass        rau0     kg/m**3
566      !!              in situ volumic mass          rho      kg/m**3
567      !!              in situ density anomalie      prd      no units
568      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
569      !!          t = 40 deg celcius, s=40 psu
570      !!      nn_eos = 1 : linear equation of state function of temperature only
571      !!              prd(t) = 0.0285 - ralpha * t
572      !!      nn_eos = 2 : linear equation of state function of temperature and
573      !!               salinity
574      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
575      !!      Note that no boundary condition problem occurs in this routine
576      !!      as (ptem,psal) are defined over the whole domain.
577      !!
578      !! ** Comments on Adjoint Routine :
579      !!      Care has been taken to avoid division by zero when computing
580      !!      the inverse of the square root of salinity at masked salinity
581      !!      points.
582      !!
583      !! ** Action  :
584      !!
585      !! References :
586      !!
587      !! History :
588      !!    8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget)  - eostan.F
589      !!    9.0 ! 07-07 (K. Mogensen) Initial version based on eostan.F
590      !!        ! 08-07 (A. Vidard) bug fix in computation of prd_tl if neos=1
591      !!-----------------------------------------------------------------------
592      !! * Modules used
593      !! * Arguments
594      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts, &  ! 1 : potential temperature  [Celcius]
595      !                                                             ! 2 : salinity               [psu]
596      &                                                     pts_tl  ! 1 : TL of potential temperature  [Celcius]
597      !                                                             ! 2 : TL of salinity               [psu]
598      REAL(wp), DIMENSION(jpi,jpj ), INTENT(  out) ::   prd_tl    ! TL of in_situ density [-]
599      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::       pdep  ! depth                  [m]
600      !
601      !! * Local declarations
602      INTEGER ::  ji, jj, jk      ! dummy loop indices
603      REAL(wp) ::   &             ! temporary scalars
604         zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw,   &
605         zb, zd, zc, zaw, za, zb1, za1, zkw, zk0,               &
606         zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl,          &
607         zr4tl, zrhoptl, zetl, zbwtl,                           &
608         zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl,           &
609         zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps,                &
610         zmask
611      REAL(wp), POINTER, DIMENSION(:,:) :: zws
612      !!----------------------------------------------------------------------
613      !
614      IF( nn_timing == 1 ) CALL timing_start('eos2d')
615      !
616      CALL wrk_alloc( jpi, jpj, zws )
617      !
618      SELECT CASE ( nn_eos )
619      !
620      CASE ( 0 )               !== Jackett and McDougall (1994) formulation ==!
621      !
622#ifdef key_sp
623         zeps = 1.e-7
624#else
625         zeps = 1.e-14
626#endif
627!CDIR NOVERRCHK
628         DO jj = 1, jpjm1
629!CDIR NOVERRCHK
630            DO ji = 1, fs_jpim1   ! vector opt.
631               zws(ji,jj) = SQRT( ABS( pts(ji,jj, jp_sal) ) )
632            END DO
633         END DO
634         DO jj = 1, jpjm1
635            DO ji = 1, fs_jpim1   ! vector opt.
636
637               zmask = tmask(ji,jj,1)      ! land/sea bottom mask = surf. mask
638
639               zt = pts (ji,jj, jp_tem)            ! interpolated T
640               zs = pts (ji,jj, jp_sal)            ! interpolated S
641               zsr= zws(ji,jj)            ! square root of interpolated S
642               zh = pdep(ji,jj)            ! depth at the partial step level
643                  ! compute volumic mass pure water at atm pressure
644                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp)*zt   &
645                     & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
646                  ! seawater volumic mass atm pressure
647                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   &
648                     & -4.0899e-3_wp ) *zt+0.824493_wp
649                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp
650                  zr4= 4.8314e-4_wp
651
652                  ! potential volumic mass (reference to the surface)
653                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
654
655                  ! add the compression terms
656                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
657                  zbw= (  1.296821e-6*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
658                  zb = zbw + ze * zs
659
660                  zd = -2.042967e-2_wp
661                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
662                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
663                  za = ( zd*zsr + zc ) *zs + zaw
664
665                  zb1=   (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp
666                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)*zt-65.00517_wp ) *zt+1044.077_wp
667                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
668                  zk0= ( zb1*zsr + za1 )*zs + zkw
669
670                  ! Tangent linear part
671
672                  zttl = pts_tl(ji,jj, jp_tem)
673                  zstl = pts_tl(ji,jj, jp_sal)
674
675                  zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
676                     &   * tmask(ji,jj,1)                 * zstl
677
678                  zr1tl= ( ( ( (  5.*6.536332e-9_wp   * zt &
679                     &           -4.*1.120083e-6_wp ) * zt &
680                     &           +3.*1.001685e-4_wp ) * zt &
681                     &           -2.*9.095290e-3_wp ) * zt &
682                     &           +   6.793952e-2_wp        ) * zttl
683
684                  zr2tl= ( ( (    4.*5.3875e-9_wp   * zt &
685                     &           -3.*8.2467e-7_wp ) * zt &
686                     &           +2.*7.6438e-5_wp ) * zt &
687                     &           -   4.0899e-3_wp        ) * zttl
688
689                  zr3tl= (       -2.*1.6546e-6_wp   * zt &
690                     &           +   1.0227e-4_wp        ) * zttl
691
692                  zrhoptl=                                  zr1tl &
693                     &     + zs                           * zr2tl &
694                     &     + zsr * zs                     * zr3tl &
695                     &     + zr3 * zs                     * zsrtl &
696                     &     + (  2. * zr4 * zs + zr2              &
697                     &        + zr3 * zsr           )     * zstl
698
699                  zetl = (       -2.*3.508914e-8_wp   * zt &
700                     &           -   1.248266e-8_wp        ) * zttl
701
702                  zbwtl= (        2.*1.296821e-6_wp   * zt &
703                     &           -   5.782165e-9_wp        ) * zttl
704
705                  zbtl =                                    zbwtl &
706                     &    + zs                            * zetl  &
707                  &       + ze                            * zstl
708
709                  zctl = (       -2.*7.267926e-5_wp   * zt &
710                     &           +   2.598241e-3_wp        ) * zttl
711
712                  zawtl= ( (      3.*5.939910e-6_wp   * zt &
713                     &           +2.*2.512549e-3_wp ) * zt &
714                     &           -   0.1028859_wp          ) * zttl
715
716                  zatl =                                    zawtl &
717                  &      + zd * zs                        * zsrtl &
718                  &      + zs                             * zctl  &
719                  &      + ( zd * zsr + zc )              * zstl
720
721                  zb1tl= (       -2.*0.1909078_wp     * zt &
722                     &           +   7.390729_wp           ) * zttl
723
724                  za1tl= ( (      3.*2.326469e-3_wp   * zt &
725                     &           +2.*1.553190_wp    ) * zt &
726                     &           -   65.00517_wp           ) * zttl
727
728                  zkwtl= ( ( (   -4.*1.361629e-4_wp   * zt &
729                     &           -3.*1.852732e-2_wp ) * zt &
730                     &           -2.*30.41638_wp    ) * zt &
731                     &           +   2098.925_wp           ) * zttl
732
733                  zk0tl=                                    zkwtl &
734                     &  + zb1 * zs                        * zsrtl &
735                     &  + zs  * zsr                       * zb1tl &
736                     &  + zs                              * za1tl &
737                     &  + ( zb1 * zsr + za1 )             * zstl
738
739                  ! Masked in situ density anomaly
740
741                  zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
742                  zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
743
744                  prd_tl(ji,jj) = tmask(ji,jj,1)  * zrdc2 *                &
745                     &            (                              zrhoptl   &
746                     &              - zrdc2 * zh * zrdc1**2 * zrhop        &
747                     &                * (                        zk0tl     &
748                     &                    - zh * (               zatl      &
749                     &                             - zh        * zbtl ) ) )&
750                     &               / rau0
751
752
753            END DO
754         END DO
755         !
756      CASE ( 1 )                !==  Linear formulation = F( temperature )  ==!
757         DO jj = 1, jpjm1
758            DO ji = 1, fs_jpim1   ! vector opt.
759               prd_tl(ji,jj) = ( - rn_alpha * pts_tl(ji,jj,jp_tem) )     * tmask(ji,jj,1)
760            END DO
761         END DO
762         !
763      CASE ( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
764         DO jj = 1, jpjm1
765            DO ji = 1, fs_jpim1   ! vector opt.
766               prd_tl (ji,jj) = ( rn_beta * pts_tl(ji,jj, jp_sal) - rn_alpha * pts_tl(ji,jj, jp_tem) ) * tmask(ji,jj,1)
767            END DO
768         END DO
769         !
770      END SELECT
771      CALL wrk_dealloc( jpi, jpj, zws )
772      !
773   END SUBROUTINE eos_insitu_2d_tan
774
775   SUBROUTINE eos_insitu_adj(pts, pts_ad, prd_ad)
776      !!-----------------------------------------------------------------------
777      !!
778      !!        ***  ROUTINE eos_insitu_tan : Adjoint OF ROUTINE eos_insitu ***
779      !!
780      !! ** Purpose of direct routine   : Compute the in situ density
781      !!       (ratio rho/rau0) from potential temperature and salinity
782      !!       using an equation of state defined through the namelist
783      !!       parameter nneos.
784      !!
785      !! ** Method of direct routine    : 3 cases:
786      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
787      !!         the in situ density is computed directly as a function of
788      !!         potential temperature relative to the surface (the opa t
789      !!         variable), salt and pressure (assuming no pressure variation
790      !!         along geopotential surfaces, i.e. the pressure p in decibars
791      !!         is approximated by the depth in meters.
792      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
793      !!         with pressure                      p        decibars
794      !!              potential temperature         t        deg celsius
795      !!              salinity                      s        psu
796      !!              reference volumic mass        rau0     kg/m**3
797      !!              in situ volumic mass          rho      kg/m**3
798      !!              in situ density anomalie      prd      no units
799      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
800      !!          t = 40 deg celcius, s=40 psu
801      !!      nn_eos = 1 : linear equation of state function of temperature only
802      !!              prd(t) = 0.0285 - rn_alpha * t
803      !!      nn_eos = 2 : linear equation of state function of temperature and
804      !!               salinity
805      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
806      !!      Note that no boundary condition problem occurs in this routine
807      !!      as (ptem,psal) are defined over the whole domain.
808      !!
809      !! ** Comments on Adjoint Routine :
810      !!      Care has been taken to avoid division by zero when computing
811      !!      the inverse of the square root of salinity at masked salinity
812      !!      points.
813      !!
814      !! ** Action  :
815      !!
816      !! References :
817      !!
818      !! History :
819      !!    8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget)  - eostan.F
820      !!    9.0 ! 08-08 (A. Vidard) 9.0 version
821      !!-----------------------------------------------------------------------
822      !! * Modules used
823      !! * Arguments
824      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts      ! 1 : potential temperature  [Celcius]
825      !                                                         ! 2 : salinity               [psu]
826      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::  pts_ad    ! 1 : TL of potential temperature [Celsius]
827                                                                ! 2 : TL of salinity [psu]
828      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   &
829         & prd_ad                   ! TL of potential density (surface referenced)
830
831      !! * Local declarations
832      INTEGER ::  ji, jj, jk      ! dummy loop indices
833      REAL(wp) ::   &             ! temporary scalars
834         zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw,   &
835         zb, zd, zc, zaw, za, zb1, za1, zkw, zk0,               &
836         ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad,          &
837         zr4ad, zrhopad, zead, zbwad,                           &
838         zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad,           &
839         zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps,                &
840         zmask, zrau0r
841      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
842      !!----------------------------------------------------------------------
843      IF( nn_timing == 1 ) CALL timing_start('eos_adj')
844      !
845      CALL wrk_alloc( jpi, jpj, jpk, zws )
846      !
847      ! initialization of adjoint variables
848      ztad    = 0.0_wp
849      zsad    = 0.0_wp
850      zhad    = 0.0_wp
851      zsrad   = 0.0_wp
852      zr1ad   = 0.0_wp
853      zr2ad   = 0.0_wp
854      zr3ad   = 0.0_wp
855      zr4ad   = 0.0_wp
856      zrhopad = 0.0_wp
857      zead    = 0.0_wp
858      zbwad   = 0.0_wp
859      zbad    = 0.0_wp
860      zdad    = 0.0_wp
861      zcad    = 0.0_wp
862      zawad   = 0.0_wp
863      zaad    = 0.0_wp
864      zb1ad   = 0.0_wp
865      za1ad   = 0.0_wp
866      zkwad   = 0.0_wp
867      zk0ad   = 0.0_wp
868      SELECT CASE ( nn_eos )
869
870      CASE ( 0 )               !== Jackett and McDougall (1994) formulation ==!
871         zrau0r = 1.e0 / rau0
872#ifdef key_sp
873         zeps = 1.e-7
874#else
875         zeps = 1.e-14
876#endif
877!CDIR NOVERRCHK
878         zws(:,:,:) = SQRT( ABS( pts(:,:,:, jp_sal) ) )
879         DO jk = jpkm1, 1, -1
880            DO jj = jpj, 1, -1
881               DO ji = jpi, 1, -1
882                  zt = pts(ji,jj,jk, jp_tem)
883                  zs = pts(ji,jj,jk, jp_sal)
884                  zh = fsdept(ji,jj,jk)        ! depth
885                  zsr= zws(ji,jj,jk)           ! square root salinity
886                  ! compute volumic mass pure water at atm pressure
887                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp)*zt   &
888                     -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
889                  ! seawater volumic mass atm pressure
890                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   &
891                     -4.0899e-3_wp ) *zt+0.824493_wp
892                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp
893                  zr4= 4.8314e-4_wp
894
895                  ! potential volumic mass (reference to the surface)
896                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
897
898                  ! add the compression terms
899                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
900                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
901                  zb = zbw + ze * zs
902
903                  zd = -2.042967e-2_wp
904                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
905                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
906                  za = ( zd*zsr + zc ) *zs + zaw
907
908                  zb1=   (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp
909                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)*zt-65.00517_wp ) *zt+1044.077_wp
910                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
911                  zk0= ( zb1*zsr + za1 )*zs + zkw
912
913                  zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
914                  zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
915                  ! ============
916                  ! Adjoint part
917                  ! ============
918
919                  ! Masked in situ density anomaly
920
921                  zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk)    &
922                     &                                 * zrdc2 * zrau0r
923                  zk0ad   = zk0ad   - prd_ad(ji,jj,jk) * tmask(ji,jj,jk)    &
924                     &                                 * zrdc2 * zrdc2 * zh &
925                     &                                 * zrdc1**2 * zrhop   &
926                     &                                 * zrau0r
927                  zaad    = zaad    + prd_ad(ji,jj,jk) * tmask(ji,jj,jk)    &
928                     &                                 * zrdc2 * zrdc2 * zh &
929                     &                                 * zrdc1**2 * zrhop   &
930                     &                                 * zh * zrau0r
931                  zbad    = zbad    - prd_ad(ji,jj,jk) * tmask(ji,jj,jk)    &
932                     &                                 * zrdc2 * zrdc2 * zh &
933                     &                                 * zrdc1**2 * zrhop   &
934                     &                                 * zh * zh * zrau0r
935                  prd_ad(ji,jj,jk) = 0.0_wp
936
937                  zkwad = zkwad + zk0ad
938                  zsrad = zsrad + zk0ad * zb1 * zs
939                  zb1ad = zb1ad + zk0ad * zs  * zsr
940                  za1ad = za1ad + zk0ad * zs
941                  zsad  = zsad  + zk0ad * ( zb1 * zsr + za1 )
942                  zk0ad = 0.0_wp
943
944                  ztad  = ztad + zkwad * ( ( (-4.*1.361629e-4_wp   * zt &
945                     &                        -3.*1.852732e-2_wp ) * zt &
946                     &                        -2.*30.41638_wp    ) * zt &
947                     &                        +   2098.925_wp           )
948                  zkwad = 0.0_wp
949
950                  ztad  = ztad + za1ad * ( ( 3.*2.326469e-3_wp   * zt &
951                     &                      +2.*1.553190_wp    ) * zt &
952                     &                      -   65.00517_wp           )
953                  za1ad = 0.0_wp
954
955                  ztad  = ztad + zb1ad * (-2.*0.1909078_wp     * zt &
956                     &                    +   7.390729_wp           )
957                  zb1ad = 0.0_wp
958
959                  zawad = zawad + zaad
960                  zsrad = zsrad + zaad *   zd * zs
961                  zcad  = zcad  + zaad *   zs
962                  zsad  = zsad  + zaad * ( zd * zsr + zc )
963                  zaad  = 0.0_wp
964
965                  ztad  = ztad + zawad * ( ( 3.*5.939910e-6_wp   * zt &
966                     &                      +2.*2.512549e-3_wp ) * zt &
967                     &                      -   0.1028859_wp          )
968                  zawad = 0.0_wp
969
970                  ztad  = ztad + zcad * (-2.*7.267926e-5_wp   * zt &
971                     &                   +   2.598241e-3_wp        )
972                  zcad  = 0.0_wp
973
974                  zbwad = zbwad + zbad
975                  zead  = zead  + zbad * zs
976                  zsad  = zsad  + zbad * ze
977                  zbad  = 0.0_wp
978
979                  ztad  = ztad + zbwad *  ( 2.*1.296821e-6_wp   * zt &
980                     &                     -   5.782165e-9_wp        )
981                  zbwad = 0.0_wp
982
983                  ztad  = ztad + zead * (-2.*3.508914e-8_wp   * zt &
984                     &                   -   1.248266e-8_wp        )
985                  zead =  0.0_wp
986
987                  zr1ad   = zr1ad + zrhopad
988                  zr2ad   = zr2ad + zrhopad * zs
989                  zr3ad   = zr3ad + zrhopad * zsr * zs
990                  zsrad   = zsrad + zrhopad * zr3 * zs
991                  zsad    = zsad  + zrhopad * ( 2. * zr4 * zs + zr2  &
992                     &                        + zr3 * zsr           )
993                  zrhopad = 0.0_wp
994
995                  ztad  = ztad + zr3ad * (-2.*1.6546e-6_wp   * zt &
996                     &                    +   1.0227e-4_wp        )
997                  zr3ad = 0.0_wp
998
999                  ztad  = ztad + zr2ad * ( ( ( 4.*5.3875e-9_wp   * zt &
1000                     &                        -3.*8.2467e-7_wp ) * zt &
1001                     &                        +2.*7.6438e-5_wp ) * zt &
1002                     &                        -   4.0899e-3_wp        )
1003                  zr2ad = 0.0_wp
1004
1005                  ztad  = ztad + zr1ad * ( ( ( ( 5.*6.536332e-9_wp   * zt &
1006                     &                          -4.*1.120083e-6_wp ) * zt &
1007                     &                          +3.*1.001685e-4_wp ) * zt &
1008                     &                          -2.*9.095290e-3_wp ) * zt &
1009                     &                          +   6.793952e-2_wp        )
1010                  zr1ad = 0.0_wp
1011
1012                  zsad  = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1013                     &                 * tmask(ji,jj,jk)
1014                  zsrad = 0.0_wp
1015
1016                  pts_ad(ji,jj,jk, jp_sal) = pts_ad(ji,jj,jk, jp_sal) + zsad
1017                  pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk, jp_tem) + ztad
1018                  ztad = 0.0_wp
1019                  zsad = 0.0_wp
1020               END DO
1021            END DO
1022         END DO
1023         !
1024      CASE ( 1 )               !== Linear formulation function of temperature only ==!
1025         DO jk = jpkm1, 1, -1
1026            pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1027            prd_ad(:,:,jk) = 0.0_wp
1028         END DO
1029         !
1030      CASE ( 2 )               !== Linear formulation function of temperature and salinity ==!
1031         DO jk = jpkm1, 1, -1
1032            pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1033            pts_ad(:,:,jk,jp_sal) = pts_ad(:,:,jk,jp_sal) + rn_beta * prd_ad( :,:,jk) * tmask(:,:,jk)
1034            prd_ad( :,:,jk) = 0.0_wp
1035         END DO
1036         !
1037      END SELECT
1038      !
1039      CALL wrk_dealloc( jpi, jpj, jpk, zws )
1040      !
1041      IF( nn_timing == 1 ) CALL timing_stop('eos_adj')
1042      !
1043   END SUBROUTINE eos_insitu_adj
1044
1045   SUBROUTINE eos_insitu_pot_adj ( pts, pts_ad, prd_ad, prhop_ad )
1046      !!----------------------------------------------------------------------
1047      !!                  ***  ROUTINE eos_insitu_pot_adj  ***
1048      !!
1049      !! ** Purpose or the direct routine:
1050      !!      Compute the in situ density (ratio rho/rau0) and the
1051      !!      potential volumic mass (Kg/m3) from potential temperature and
1052      !!      salinity fields using an equation of state defined through the
1053      !!     namelist parameter nn_eos.
1054      !!
1055      !! ** Method  :
1056      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
1057      !!         the in situ density is computed directly as a function of
1058      !!         potential temperature relative to the surface (the opa t
1059      !!         variable), salt and pressure (assuming no pressure variation
1060      !!         along geopotential surfaces, i.e. the pressure p in decibars
1061      !!         is approximated by the depth in meters.
1062      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
1063      !!              rhop(t,s)  = rho(t,s,0)
1064      !!         with pressure                      p        decibars
1065      !!              potential temperature         t        deg celsius
1066      !!              salinity                      s        psu
1067      !!              reference volumic mass        rau0     kg/m**3
1068      !!              in situ volumic mass          rho      kg/m**3
1069      !!              in situ density anomalie      prd      no units
1070      !!
1071      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
1072      !!          t = 40 deg celcius, s=40 psu
1073      !!
1074      !!      neos = 1 : linear equation of state function of temperature only
1075      !!              prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - ralpha * t
1076      !!              rhop(t,s)  = rho(t,s)
1077      !!
1078      !!      nn_eos = 2 : linear equation of state function of temperature and
1079      !!               salinity
1080      !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0
1081      !!                       = rn_beta * s - rn_alpha * tn - 1.
1082      !!              rhop(t,s)  = rho(t,s)
1083      !!      Note that no boundary condition problem occurs in this routine
1084      !!      as (tn,sn) or (ta,sa) are defined over the whole domain.
1085      !!
1086      !! ** Action  : - prd  , the in situ density (no units)
1087      !!              - prhop, the potential volumic mass (Kg/m3)
1088      !!
1089      !! References :
1090      !!      Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
1091      !!      Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
1092      !!
1093      !! History of the adjoint routine:
1094      !!   9.0  !  08-06  (A. Vidard) Initial version
1095      !!----------------------------------------------------------------------
1096      !! * Arguments
1097
1098      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts ! 1 : potential temperature/salinity [Celcius/psu]
1099      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) ::   pts_ad ! 1 : potential temperature/salinity [Celcius/psu]
1100      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(inout) ::   prd_ad      ! TL of in_situ density [-]
1101      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(inout) ::   prhop_ad    ! TL of potential density (surface referenced)
1102      !! * Local declarations
1103      INTEGER ::  ji, jj, jk      ! dummy loop indices
1104      REAL(wp) ::   &             ! temporary scalars
1105         zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw,   &
1106         zb, zd, zc, zaw, za, zb1, za1, zkw, zk0,               &
1107         ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad,          &
1108         zr4ad, zrhopad, zead, zbwad,                           &
1109         zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad,           &
1110         zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps,                &
1111         zmask, zrau0r
1112      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
1113      !!----------------------------------------------------------------------
1114      !
1115      IF( nn_timing == 1 ) CALL timing_start('eos-p_adj')
1116      !
1117      CALL wrk_alloc( jpi, jpj, jpk, zws )
1118      !
1119      ! initialization of adjoint variables
1120      ztad = 0.0_wp
1121      zsad = 0.0_wp
1122      zhad = 0.0_wp
1123      zsrad = 0.0_wp
1124      zr1ad = 0.0_wp
1125      zr2ad = 0.0_wp
1126      zr3ad = 0.0_wp
1127      zr4ad = 0.0_wp
1128      zrhopad = 0.0_wp
1129      zead = 0.0_wp
1130      zbwad = 0.0_wp
1131      zbad = 0.0_wp
1132      zdad = 0.0_wp
1133      zcad = 0.0_wp
1134      zawad = 0.0_wp
1135      zaad = 0.0_wp
1136      zb1ad  = 0.0_wp
1137      za1ad = 0.0_wp
1138      zkwad = 0.0_wp
1139      zk0ad = 0.0_wp
1140
1141      SELECT CASE ( nn_eos )
1142
1143      CASE ( 0 )               !== Jackett and McDougall (1994) formulation ==!
1144         zrau0r = 1.e0 / rau0
1145#ifdef key_sp
1146         zeps = 1.e-7
1147#else
1148         zeps = 1.e-14
1149#endif
1150!CDIR NOVERRCHK
1151         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
1152         !
1153         DO jk =  jpkm1, 1, -1
1154            DO jj = jpj, 1, -1
1155               DO ji = jpi, 1, -1
1156                  ! direct recomputing
1157                  zt = pts(ji,jj,jk,jp_tem)
1158                  zs = pts(ji,jj,jk,jp_sal)
1159                  zh = fsdept(ji,jj,jk)        ! depth
1160                  zsr = zws(ji,jj,jk)          ! square root salinity
1161                  ! compute volumic mass pure water at atm pressure
1162                  zr1 = ( ( ( ( 6.536332e-9_wp   * zt - 1.120083e-6_wp ) * zt &
1163                     &        + 1.001685e-4_wp ) * zt - 9.095290e-3_wp ) * zt &
1164                     &        + 6.793952e-2_wp ) * zt + 999.842594_wp
1165                  ! seawater volumic mass atm pressure
1166                  zr2 = ( ( ( 5.3875e-9_wp   * zt - 8.2467e-7_wp ) * zt       &
1167                     &      + 7.6438e-5_wp ) * zt - 4.0899e-3_wp ) * zt + 0.824493_wp
1168                  zr3 = ( -1.6546e-6_wp * zt + 1.0227e-4_wp ) * zt - 5.72466e-3_wp
1169                  zr4 = 4.8314e-4_wp
1170                  ! potential volumic mass (reference to the surface)
1171                  zrhop = ( zr4 * zs + zr3*zsr + zr2 ) * zs + zr1
1172                  ! add the compression terms
1173                  ze  = ( -3.508914e-8_wp * zt - 1.248266e-8_wp ) * zt - 2.595994e-6_wp
1174                  zbw = (  1.296821e-6_wp * zt - 5.782165e-9_wp ) * zt + 1.045941e-4_wp
1175                  zb  = zbw + ze * zs
1176
1177                  zd = -2.042967e-2_wp
1178                  zc =   (-7.267926e-5_wp * zt + 2.598241e-3_wp ) * zt + 0.1571896_wp
1179                  zaw= ( ( 5.939910e-6_wp * zt + 2.512549e-3_wp ) * zt - 0.1028859_wp &
1180                     &   ) * zt - 4.721788
1181                  za = ( zd * zsr + zc ) * zs + zaw
1182
1183                  zb1=   (-0.1909078_wp   * zt + 7.390729_wp ) * zt - 55.87545_wp
1184                  za1= ( ( 2.326469e-3_wp * zt + 1.553190_wp ) * zt - 65.00517_wp &
1185                     &   ) * zt + 1044.077_wp
1186                  zkw= ( ( (-1.361629e-4_wp * zt - 1.852732e-2_wp ) * zt - 30.41638_wp &
1187                     &    ) * zt + 2098.925_wp ) * zt + 190925.6_wp
1188                  zk0= ( zb1 * zsr + za1 ) * zs + zkw
1189
1190
1191                  zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1192                  zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1193
1194                  ! ============
1195                  ! Adjoint part
1196                  ! ============
1197
1198                  ! Masked in situ density anomaly
1199
1200                  zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk)    &
1201                     &                                 * zrdc2 * zrau0r
1202                  zk0ad   = zk0ad   - prd_ad(ji,jj,jk) * tmask(ji,jj,jk)    &
1203                     &                                 * zrdc2 * zrdc2 * zh &
1204                     &                                 * zrdc1**2 * zrhop   &
1205                     &                                 * zrau0r
1206                  zaad    = zaad    + prd_ad(ji,jj,jk) * tmask(ji,jj,jk)    &
1207                     &                                 * zrdc2 * zrdc2 * zh &
1208                     &                                 * zrdc1**2 * zrhop   &
1209                     &                                 * zh * zrau0r
1210                  zbad    = zbad    - prd_ad(ji,jj,jk) * tmask(ji,jj,jk)    &
1211                     &                                 * zrdc2 * zrdc2 * zh &
1212                     &                                 * zrdc1**2 * zrhop   &
1213                     &                                 * zh * zh * zrau0r
1214                  prd_ad(ji,jj,jk) = 0.0_wp
1215
1216                  zkwad = zkwad + zk0ad
1217                  zsrad = zsrad + zk0ad * zb1 * zs
1218                  zb1ad = zb1ad + zk0ad * zs  * zsr
1219                  za1ad = za1ad + zk0ad * zs
1220                  zsad  = zsad  + zk0ad * ( zb1 * zsr + za1 )
1221                  zk0ad = 0.0_wp
1222
1223                  ztad  = ztad + zkwad * ( ( (-4.*1.361629e-4_wp   * zt &
1224                     &                        -3.*1.852732e-2_wp ) * zt &
1225                     &                        -2.*30.41638_wp    ) * zt &
1226                     &                        +   2098.925_wp           )
1227                  zkwad = 0.0_wp
1228
1229                  ztad  = ztad + za1ad * ( ( 3.*2.326469e-3_wp   * zt &
1230                     &                      +2.*1.553190_wp    ) * zt &
1231                     &                      -   65.00517_wp           )
1232                  za1ad = 0.0_wp
1233
1234                  ztad  = ztad + zb1ad * (-2.*0.1909078_wp     * zt &
1235                     &                    +   7.390729_wp           )
1236                  zb1ad = 0.0_wp
1237
1238                  zawad = zawad + zaad
1239                  zsrad = zsrad + zaad *   zd * zs
1240                  zcad  = zcad  + zaad *   zs
1241                  zsad  = zsad  + zaad * ( zd * zsr + zc )
1242                  zaad  = 0.0_wp
1243
1244                  ztad  = ztad + zawad * ( ( 3.*5.939910e-6_wp   * zt &
1245                     &                      +2.*2.512549e-3_wp ) * zt &
1246                     &                      -   0.1028859_wp          )
1247                  zawad = 0.0_wp
1248
1249                  ztad  = ztad + zcad * (-2.*7.267926e-5_wp   * zt &
1250                     &                   +   2.598241e-3_wp        )
1251                  zcad  = 0.0_wp
1252
1253
1254                  zsad  = zsad  + zbad * ze
1255                  zead  = zead  + zbad * zs
1256                  zbwad = zbwad + zbad
1257                  zbad  = 0.0_wp
1258
1259                  ztad  = ztad + zbwad *  ( 2.*1.296821e-6_wp   * zt &
1260                     &                     -   5.782165e-9_wp        )
1261                  zbwad = 0.0_wp
1262
1263                  ztad  = ztad + zead * (-2.*3.508914e-8_wp   * zt &
1264                     &                   -   1.248266e-8_wp        )
1265                  zead =  0.0_wp
1266
1267                  zrhopad = zrhopad + prhop_ad(ji,jj,jk) * tmask(ji,jj,jk)
1268                  prhop_ad(ji,jj,jk) = 0.0_wp
1269
1270                  zr1ad   = zr1ad + zrhopad
1271                  zr2ad   = zr2ad + zrhopad * zs
1272                  zr3ad   = zr3ad + zrhopad * zsr * zs
1273                  zsrad   = zsrad + zrhopad * zr3 * zs
1274                  zsad    = zsad  + zrhopad * ( 2. * zr4 * zs + zr2  &
1275                     &                        + zr3 * zsr           )
1276                  zrhopad = 0.0_wp
1277
1278                  ztad  = ztad + zr3ad * (-2.*1.6546e-6_wp   * zt &
1279                     &                    +   1.0227e-4_wp        )
1280                  zr3ad = 0.0_wp
1281
1282                  ztad  = ztad + zr2ad * ( ( ( 4.*5.3875e-9_wp   * zt &
1283                     &                        -3.*8.2467e-7_wp ) * zt &
1284                     &                        +2.*7.6438e-5_wp ) * zt &
1285                     &                        -   4.0899e-3_wp        )
1286                  zr2ad = 0.0_wp
1287
1288                  ztad  = ztad + zr1ad * ( ( ( ( 5.*6.536332e-9_wp   * zt &
1289                     &                          -4.*1.120083e-6_wp ) * zt &
1290                     &                          +3.*1.001685e-4_wp ) * zt &
1291                     &                          -2.*9.095290e-3_wp ) * zt &
1292                     &                          +   6.793952e-2_wp        )
1293                  zr1ad = 0.0_wp
1294
1295                  zsad  = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1296                     &                 * tmask(ji,jj,jk)
1297                  zsrad = 0.0_wp
1298
1299                  pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + zsad
1300                  pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + ztad
1301                  ztad = 0.0_wp
1302                  zsad = 0.0_wp
1303               END DO
1304            END DO
1305         END DO
1306         !
1307      CASE ( 1 )               !==  Linear formulation = F( temperature )  ==!
1308         DO jk = jpkm1, 1, -1
1309            prd_ad(:,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk)
1310            prhop_ad(:,:,jk) = 0.0_wp
1311            pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) -  rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1312            prd_ad(:,:,jk) = 0.0_wp
1313         END DO
1314         !
1315      CASE ( 2 )               !==  Linear formulation = F( temperature , salinity )  ==!
1316         DO jk = jpkm1, 1, -1
1317            prd_ad(  :,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk)
1318            prhop_ad(:,:,jk) = 0.0_wp
1319            pts_ad( :,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1320            pts_ad( :,:,jk,jp_sal) = pts_ad(:,:,jk,jp_sal) + rn_beta   * prd_ad(:,:,jk) * tmask(:,:,jk)
1321            prd_ad(  :,:,jk) = 0.0_wp
1322         END DO
1323         !
1324      END SELECT
1325      CALL wrk_dealloc( jpi, jpj, jpk, zws )
1326      !
1327      IF( nn_timing == 1 ) CALL timing_stop('eos-p_adj')
1328      !
1329   END SUBROUTINE eos_insitu_pot_adj
1330
1331   SUBROUTINE eos_insitu_2d_adj( pts, pdep, pts_ad, prd_ad )
1332      !!-----------------------------------------------------------------------
1333      !!
1334      !!        ***  ROUTINE eos_insitu_2d_adj : adj OF ROUTINE eos_insitu_2d ***
1335      !!
1336      !! ** Purpose of direct routine   : Compute the in situ density
1337      !!       (ratio rho/rau0) from potential temperature and salinity
1338      !!       using an equation of state defined through the namelist
1339      !!       parameter nn_eos. * 2D field case
1340      !!
1341      !! ** Method of direct routine    : 3 cases:
1342      !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state.
1343      !!         the in situ density is computed directly as a function of
1344      !!         potential temperature relative to the surface (the opa t
1345      !!         variable), salt and pressure (assuming no pressure variation
1346      !!         along geopotential surfaces, i.e. the pressure p in decibars
1347      !!         is approximated by the depth in meters.
1348      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
1349      !!         with pressure                      p        decibars
1350      !!              potential temperature         t        deg celsius
1351      !!              salinity                      s        psu
1352      !!              reference volumic mass        rau0     kg/m**3
1353      !!              in situ volumic mass          rho      kg/m**3
1354      !!              in situ density anomalie      prd      no units
1355      !!         Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
1356      !!          t = 40 deg celcius, s=40 psu
1357      !!      nn_eos = 1 : linear equation of state function of temperature only
1358      !!              prd(t) = 0.0285 - rn_alpha * t
1359      !!      nn_eos = 2 : linear equation of state function of temperature and
1360      !!               salinity
1361      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
1362      !!      Note that no boundary condition problem occurs in this routine
1363      !!      as (ptem,psal) are defined over the whole domain.
1364      !!
1365      !! ** Comments on Adjoint Routine :
1366      !!      Care has been taken to avoid division by zero when computing
1367      !!      the inverse of the square root of salinity at masked salinity
1368      !!      points.
1369      !!
1370      !! ** Action  :
1371      !!
1372      !! References :
1373      !!
1374      !! History :
1375      !!    8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget)  - eosadj.F
1376      !!    9.0 ! 08-07 (A. Vidard) first version based on eosadj
1377      !!-----------------------------------------------------------------------
1378      !! * Modules used
1379      !! * Arguments
1380      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius]
1381      !                                                            ! 2 : salinity               [psu]
1382       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) ::  pts_ad ! 1 : TL of potential temperature  [Celcius]
1383      !                                                            ! 2 : TL of salinity               [psu]
1384      REAL(wp), DIMENSION(jpi,jpj ), INTENT(  inout) ::    prd_ad  ! TL of in_situ density [-]
1385      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::    pdep  ! depth                  [m]
1386      !
1387
1388      INTEGER ::  ji, jj          ! dummy loop indices
1389      REAL(wp) ::   &             ! temporary scalars
1390         zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw,   &
1391         zb, zd, zc, zaw, za, zb1, za1, zkw, zk0,               &
1392         ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad,          &
1393         zr4ad, zrhopad, zead, zbwad,                           &
1394         zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad,           &
1395         zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps,                &
1396         zmask
1397      REAL(wp), POINTER, DIMENSION(:,:) :: zws
1398      !!----------------------------------------------------------------------
1399      !
1400      IF( nn_timing == 1 ) CALL timing_start('eos2d_adj')
1401      !
1402      CALL wrk_alloc( jpi, jpj, zws )
1403      !
1404      ! initialization of adjoint variables
1405      ztad    = 0.0_wp
1406      zsad    = 0.0_wp
1407      zhad    = 0.0_wp
1408      zsrad   = 0.0_wp
1409      zr1ad   = 0.0_wp
1410      zr2ad   = 0.0_wp
1411      zr3ad   = 0.0_wp
1412      zr4ad   = 0.0_wp
1413      zrhopad = 0.0_wp
1414      zead    = 0.0_wp
1415      zbwad   = 0.0_wp
1416      zbad    = 0.0_wp
1417      zdad    = 0.0_wp
1418      zcad    = 0.0_wp
1419      zawad   = 0.0_wp
1420      zaad    = 0.0_wp
1421      zb1ad   = 0.0_wp
1422      za1ad   = 0.0_wp
1423      zkwad   = 0.0_wp
1424      zk0ad   = 0.0_wp
1425      SELECT CASE ( nn_eos )
1426      CASE ( 0 )               !== Jackett and McDougall (1994) formulation ==!
1427
1428#ifdef key_sp
1429         zeps = 1.e-7
1430#else
1431         zeps = 1.e-14
1432#endif
1433!CDIR NOVERRCHK
1434         DO jj = jpjm1, 1, -1
1435!CDIR NOVERRCHK
1436            DO ji = fs_jpim1, 1, -1   ! vector opt.
1437               zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
1438            END DO
1439         END DO
1440         !
1441         DO jj = jpjm1, 1, -1
1442            DO ji = fs_jpim1, 1, -1   ! vector opt.
1443               zmask = tmask(ji,jj,1)      ! land/sea bottom mask = surf. mask
1444               zt = pts (ji,jj,jp_tem)           ! interpolated T
1445               zs = pts (ji,jj,jp_sal)           ! interpolated S
1446               zsr= zws(ji,jj)             ! square root of interpolated S
1447               zh = pdep(ji,jj)            ! depth at the partial step level
1448                  ! compute volumic mass pure water at atm pressure
1449                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp)*zt   &
1450                     &  -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
1451                  ! seawater volumic mass atm pressure
1452                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   &
1453                     &  -4.0899e-3_wp ) *zt+0.824493_wp
1454                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp
1455                  zr4= 4.8314e-4_wp
1456
1457                  ! potential volumic mass (reference to the surface)
1458                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1
1459
1460                  ! add the compression terms
1461                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
1462                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
1463                  zb = zbw + ze * zs
1464
1465                  zd = -2.042967e-2_wp
1466                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
1467                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
1468                  za = ( zd*zsr + zc ) *zs + zaw
1469
1470                  zb1=   (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp
1471                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)*zt-65.00517_wp ) *zt+1044.077_wp
1472                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
1473                  zk0= ( zb1*zsr + za1 )*zs + zkw
1474
1475                  zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1476                  zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1477                  ! ============
1478                  ! Adjoint part
1479                  ! ============
1480
1481                  ! Masked in situ density anomaly
1482
1483                  zrhopad = zrhopad + prd_ad(ji,jj) * tmask(ji,jj,1)     &
1484                     &                              * zrdc2 / rau0
1485                  zk0ad   = zk0ad   - prd_ad(ji,jj) * tmask(ji,jj,1)     &
1486                     &                              * zrdc2 * zrdc2 * zh &
1487                     &                              * zrdc1**2 * zrhop   &
1488                     &                              / rau0
1489                  zaad    = zaad    + prd_ad(ji,jj) * tmask(ji,jj,1)     &
1490                     &                              * zrdc2 * zrdc2 * zh &
1491                     &                              * zrdc1**2 * zrhop   &
1492                     &                              * zh / rau0
1493                  zbad    = zbad    - prd_ad(ji,jj) * tmask(ji,jj,1)     &
1494                     &                              * zrdc2 * zrdc2 * zh &
1495                     &                              * zrdc1**2 * zrhop   &
1496                     &                              * zh * zh / rau0
1497                  prd_ad(ji,jj) = 0.0_wp
1498
1499                  zkwad = zkwad + zk0ad
1500                  zsrad = zsrad + zk0ad * zb1 * zs
1501                  zb1ad = zb1ad + zk0ad * zs  * zsr
1502                  za1ad = za1ad + zk0ad * zs
1503                  zsad  = zsad  + zk0ad * ( zb1 * zsr + za1 )
1504                  zk0ad = 0.0_wp
1505
1506                  ztad  = ztad + zkwad * ( ( (-4.*1.361629e-4_wp   * zt &
1507                     &                        -3.*1.852732e-2_wp ) * zt &
1508                     &                        -2.*30.41638_wp    ) * zt &
1509                     &                        +   2098.925_wp           )
1510                  zkwad = 0.0_wp
1511
1512                  ztad  = ztad + za1ad * ( ( 3.*2.326469e-3_wp   * zt &
1513                     &                      +2.*1.553190_wp    ) * zt &
1514                     &                      -   65.00517_wp           )
1515                  za1ad = 0.0_wp
1516
1517                  ztad  = ztad + zb1ad * (-2.*0.1909078_wp     * zt &
1518                     &                    +   7.390729_wp           )
1519                  zb1ad = 0.0_wp
1520
1521                  zawad = zawad + zaad
1522                  zsrad = zsrad + zaad *   zd * zs
1523                  zcad  = zcad  + zaad *   zs
1524                  zsad  = zsad  + zaad * ( zd * zsr + zc )
1525                  zaad  = 0.0_wp
1526
1527                  ztad  = ztad + zawad * ( ( 3.*5.939910e-6_wp   * zt &
1528                     &                      +2.*2.512549e-3_wp ) * zt &
1529                     &                      -   0.1028859_wp          )
1530                  zawad = 0.0_wp
1531
1532                  ztad  = ztad + zcad * (-2.*7.267926e-5_wp   * zt &
1533                     &                   +   2.598241e-3_wp        )
1534                  zcad  = 0.0_wp
1535
1536                  zbwad = zbwad + zbad
1537                  zead  = zead  + zbad * zs
1538                  zsad  = zsad  + zbad * ze
1539                  zbad  = 0.0_wp
1540
1541                  ztad  = ztad + zbwad *  ( 2.*1.296821e-6_wp   * zt &
1542                     &                     -   5.782165e-9_wp        )
1543                  zbwad = 0.0_wp
1544
1545                  ztad  = ztad + zead * (-2.*3.508914e-8_wp   * zt &
1546                     &                   -   1.248266e-8_wp        )
1547                  zead =  0.0_wp
1548
1549                  zr1ad   = zr1ad + zrhopad
1550                  zr2ad   = zr2ad + zrhopad * zs
1551                  zr3ad   = zr3ad + zrhopad * zsr * zs
1552                  zsrad   = zsrad + zrhopad * zr3 * zs
1553                  zsad    = zsad  + zrhopad * ( 2. * zr4 * zs + zr2  &
1554                     &                        + zr3 * zsr           )
1555                  zrhopad = 0.0_wp
1556
1557                  ztad  = ztad + zr3ad * (-2.*1.6546e-6_wp   * zt &
1558                     &                    +   1.0227e-4_wp        )
1559                  zr3ad = 0.0_wp
1560
1561                  ztad  = ztad + zr2ad * ( ( ( 4.*5.3875e-9_wp   * zt &
1562                     &                        -3.*8.2467e-7_wp ) * zt &
1563                     &                        +2.*7.6438e-5_wp ) * zt &
1564                     &                        -   4.0899e-3_wp        )
1565                  zr2ad = 0.0_wp
1566
1567                  ztad  = ztad + zr1ad * ( ( ( ( 5.*6.536332e-9_wp   * zt &
1568                     &                          -4.*1.120083e-6_wp ) * zt &
1569                     &                          +3.*1.001685e-4_wp ) * zt &
1570                     &                          -2.*9.095290e-3_wp ) * zt &
1571                     &                          +   6.793952e-2_wp        )
1572                  zr1ad = 0.0_wp
1573
1574                  zsad  = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1575                     &                 * tmask(ji,jj, 1)
1576                  zsrad = 0.0_wp
1577
1578                  pts_ad(ji,jj,jp_sal) = pts_ad(ji,jj,jp_sal) + zsad
1579                  pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) + ztad
1580                  ztad = 0.0_wp
1581                  zsad = 0.0_wp
1582            END DO
1583         END DO
1584         !
1585      CASE ( 1 )                !==  Linear formulation = F( temperature )  ==!
1586         DO jj = jpjm1, 1, -1
1587            DO ji = fs_jpim1, 1, -1   ! vector opt.
1588               pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1)
1589               prd_ad(ji,jj) = 0.0_wp
1590            END DO
1591         END DO
1592         !
1593      CASE ( 2 )               !==  Linear formulation = F( temperature , salinity )  ==!
1594         DO jj = jpjm1, 1, -1
1595            DO ji = fs_jpim1, 1, -1   ! vector opt.
1596               pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1)
1597               pts_ad(ji,jj,jp_sal) = pts_ad(ji,jj,jp_sal) + prd_ad(ji,jj) * rn_beta  * tmask(ji,jj,1)
1598               prd_ad (ji,jj) = 0.0_wp
1599            END DO
1600         END DO
1601         !
1602      END SELECT
1603      !
1604      CALL wrk_dealloc( jpi, jpj, zws )
1605      !
1606      IF( nn_timing == 1 ) CALL timing_stop('eos2d_adj')
1607      !
1608   END SUBROUTINE eos_insitu_2d_adj
1609
1610   SUBROUTINE eos_bn2_tan ( pts, pts_tl, pn2_tl )
1611      !!----------------------------------------------------------------------
1612      !!                  ***  ROUTINE eos_bn2_tan  ***
1613      !!
1614      !! ** Purpose of the direct routine:   Compute the local
1615      !!      Brunt-Vaisala frequency at the time-step of the input arguments
1616      !!
1617      !! ** Method of the direct routine:
1618      !!       * nn_eos = 0  : UNESCO sea water properties
1619      !!         The brunt-vaisala frequency is computed using the polynomial
1620      !!      polynomial expression of McDougall (1987):
1621      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
1622      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
1623      !!      computed and used in zdfddm module :
1624      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
1625      !!       * nn_eos = 1  : linear equation of state (temperature only)
1626      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w
1627      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
1628      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
1629      !!      The use of potential density to compute N^2 introduces e r r o r
1630      !!      in the sign of N^2 at great depths. We recommand the use of
1631      !!      nn_eos = 0, except for academical studies.
1632      !!        Macro-tasked on horizontal slab (jk-loop)
1633      !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr
1634      !!      and is never used at this level.
1635      !!
1636      !! ** Action  : - pn2 : the brunt-vaisala frequency
1637      !!
1638      !! References :
1639      !!      McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
1640      !!
1641      !! History:
1642      !!        !  08-07  (A. Vidard) First version
1643      !!----------------------------------------------------------------------
1644      !! * Arguments
1645      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts, &   ! 1 : potential temperature  [Celcius]
1646      !                                                                  ! 2 : salinity               [psu]
1647      &                                                         pts_tl   ! 1 : TL of potential temperature [Celsius]
1648                                                                         ! 2 : TL of salinity [psu]
1649      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
1650         & pn2_tl                   ! TL of potential density (surface referenced)
1651      !! * Local declarations
1652      INTEGER  ::   ji, jj, jk   ! dummy loop indices
1653      REAL(wp) ::             &
1654         zgde3w, zt, zs, zh,  &  ! temporary scalars
1655         zalbet, zbeta           !    "         "
1656      REAL(wp) ::             &
1657         zttl, zstl,          &  ! temporary scalars
1658         zalbettl, zbetatl       !    "         "
1659#if defined key_zdfddm
1660      REAL(wp) ::   zds, zdstl   ! temporary scalars
1661#endif
1662
1663      ! pn2_tl : interior points only (2=< jk =< jpkm1 )
1664      ! --------------------------
1665      SELECT CASE ( nn_eos )
1666
1667      CASE ( 0 )               !== Jackett and McDougall (1994) formulation ==!
1668         DO jk = 2, jpkm1
1669            DO jj = 1, jpj
1670               DO ji = 1, jpi
1671                  zgde3w = grav / fse3w(ji,jj,jk)
1672                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )          ! potential temperature at w-point
1673                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0   ! salinity anomaly (s-35) at w-point
1674                  zh = fsdepw(ji,jj,jk)                                     ! depth in meters  at w-point
1675
1676                  zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &   ! ratio alpha/beta
1677                     &                               - 0.203814e-03_wp ) * zt   &
1678                     &                               + 0.170907e-01_wp ) * zt   &
1679                     &   + 0.665157e-01_wp                                      &
1680                     &   +     ( - 0.678662e-05_wp * zs                         &
1681                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
1682                     &   +   ( ( - 0.302285e-13_wp * zh                         &
1683                     &           - 0.251520e-11_wp * zs                         &
1684                     &           + 0.512857e-12_wp * zt * zt           ) * zh   &
1685                     &           - 0.164759e-06_wp * zs                         &
1686                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
1687                     &                               + 0.380374e-04_wp ) * zh
1688
1689                  zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta
1690                     &                            - 0.301985e-05_wp ) * zt      &
1691                     &   + 0.785567e-03_wp                                      &
1692                     &   + (     0.515032e-08_wp * zs                           &
1693                     &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     &
1694                     &   +(  (   0.121551e-17_wp * zh                           &
1695                     &         - 0.602281e-15_wp * zs                           &
1696                     &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     &
1697                     &                             + 0.408195e-10_wp   * zs     &
1698                     &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     &
1699                     &                             - 0.121555e-07_wp ) * zh
1700
1701                  !! tangent part
1702                  zttl = 0.5 * ( pts_tl(ji,jj,jk,jp_tem) + pts_tl(ji,jj,jk-1,jp_tem) )  ! potential temperature at w-point
1703                  zstl = 0.5 * ( pts_tl(ji,jj,jk,jp_sal) + pts_tl(ji,jj,jk-1,jp_sal) )  ! salinity anomaly at w-point
1704                  zalbettl = (    (   ( -4.*0.255019e-07_wp   * zt                    & ! ratio alpha/beta
1705                     &              +3.*0.298357e-05_wp ) * zt                        &
1706                     &              -2.*0.203814e-03_wp ) * zt                        &
1707                     &      +           0.170907e-01_wp                               &
1708                     &      -           0.846960e-04_wp   * zs                        &
1709                     &      - (         0.933746e-06_wp                               &
1710                     &          - (  2.*0.791325e-08_wp                               &
1711                     &              +2.*0.512857e-12_wp   * zh ) * zt ) * zh ) * zttl &
1712                     & + (  -        2.*0.678662e-05_wp   * zs                        &
1713                     &      -           0.846960e-04_wp   * zt                        &
1714                     &      +           0.378110e-02_wp                               &
1715                     &      + (     -   0.164759e-06_wp                               &
1716                     &              -   0.251520e-11_wp   * zh ) * zh         ) * zstl
1717
1718                  zbetatl = (    (     -3.*0.415613e-09_wp   * zt               &
1719                     &              +2.*0.555579e-07_wp ) * zt                  &
1720                     &      -           0.301985e-05_wp                         &
1721                     &      +           0.788212e-08_wp   * zs                  &
1722                     &      + (     -2.*0.213127e-11_wp   * zt                  &
1723                     &              -   0.175379e-14_wp   * zh                  &
1724                     &              +   0.192867e-09_wp         ) * zh ) * zttl &
1725                     & + (           2.*0.515032e-08_wp   * zs                  &
1726                     &      +           0.788212e-08_wp   * zt                  &
1727                     &      -           0.356603e-06_wp                         &
1728                     &      + (     -   0.602281e-15_wp   * zh                  &
1729                     &              +   0.408195e-10_wp         ) * zh ) * zstl
1730
1731                     pn2_tl(ji,jj,jk) = zgde3w * tmask(ji,jj,jk) * (                                   &
1732                  &       zbeta   * (  zalbet                                        &
1733                  &                  * ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) )   &
1734                  &                 +  zalbettl                                      &
1735                  &                 * ( pts  (ji,jj,jk-1,jp_tem) - pts  (ji,jj,jk,jp_tem) )      &
1736                  &                 -  ( pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal) ) ) &
1737                  &     + zbetatl * (  zalbet                                        &
1738                  &                  * ( pts  (ji,jj,jk-1,jp_tem) - pts  (ji,jj,jk,jp_tem) )     &
1739                  &                 -  ( pts  (ji,jj,jk-1,jp_sal) - pts  (ji,jj,jk,jp_sal) ) ) )
1740#if defined key_zdfddm
1741                     zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1742                     zdstl = ( pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal) )
1743                     IF ( ABS( zds) <= 1.e-20 ) THEN
1744                        zds = 1.e-20
1745                        rrau_tl(ji,jj,jk) = zalbettl * &
1746                             &    ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds       &
1747                             &              + zalbet *   &
1748                             &  ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds )
1749                     ELSE
1750                     rrau_tl(ji,jj,jk) = zalbettl * &
1751                        &    ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds       &
1752                        &              + zalbet *   &
1753                        &  ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds &
1754                        &  - ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) * zdstl/ zds**2 )
1755                  ENDIF
1756#endif
1757               END DO
1758            END DO
1759         END DO
1760         !
1761      CASE ( 1 )                !==  Linear formulation = F( temperature )  ==!
1762         DO jk = 2, jpkm1
1763            pn2_tl(:,:,jk) = rn_alpha * ( pts_tl(:,:,jk-1,jp_tem) - pts_tl(:,:,jk,jp_tem) ) * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1764         END DO
1765         !
1766      CASE ( 2 )               !==  Linear formulation = F( temperature , salinity )  ==!
1767         DO jk = 2, jpkm1
1768            pn2_tl(:,:,jk) = ( rn_alpha * ( pts_tl(:,:,jk-1,jp_tem) - pts_tl(:,:,jk,jp_tem) )    &
1769                             & - rn_beta  * ( pts_tl(:,:,jk-1,jp_sal) - pts_tl(:,:,jk,jp_sal) )  ) &
1770                             & * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1771         END DO
1772#if defined key_zdfddm
1773         DO jk = 2, jpkm1
1774            DO jj = 1, jpj
1775               DO ji = 1, jpi
1776                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1777                  zdstl = pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal)
1778                  IF ( ABS( zds) <= 1.e-20 ) THEN
1779                     zds = 1.e-20
1780                     rrau_tl(ji,jj,jk) = ralpbet * &
1781                     & ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds )
1782                  ELSE
1783                     rrau_tl(ji,jj,jk) = ralpbet * &
1784                     & ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds &
1785                     & - ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) * zdstl / zds**2 )
1786                  ENDIF
1787                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1788               END DO
1789            END DO
1790         END DO
1791#endif
1792      END SELECT
1793   END SUBROUTINE eos_bn2_tan
1794
1795   SUBROUTINE eos_bn2_adj ( pts, pts_ad, pn2_ad )
1796      !!----------------------------------------------------------------------
1797      !!                  ***  ROUTINE eos_bn2_adj  ***
1798      !!
1799      !! ** Purpose of the direct routine:   Compute the local
1800      !!      Brunt-Vaisala frequency at the time-step of the input arguments
1801      !!
1802      !! ** Method of the direct routine:
1803      !!       * nn_eos = 0  : UNESCO sea water properties
1804      !!         The brunt-vaisala frequency is computed using the polynomial
1805      !!      polynomial expression of McDougall (1987):
1806      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
1807      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
1808      !!      computed and used in zdfddm module :
1809      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
1810      !!       * nn_eos = 1  : linear equation of state (temperature only)
1811      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w
1812      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
1813      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
1814      !!      The use of potential density to compute N^2 introduces e r r o r
1815      !!      in the sign of N^2 at great depths. We recommand the use of
1816      !!      nn_eos = 0, except for academical studies.
1817      !!        Macro-tasked on horizontal slab (jk-loop)
1818      !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr
1819      !!      and is never used at this level.
1820      !!
1821      !! ** Action  : - pn2 : the brunt-vaisala frequency
1822      !!
1823      !! References :
1824      !!      McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
1825      !!
1826      !! History:
1827      !!        !  08-07  (A. Vidard) First version
1828      !!----------------------------------------------------------------------
1829      !! * Arguments
1830      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts       ! 1 : potential temperature  [Celcius]
1831      !                                                                  ! 2 : salinity               [psu]
1832      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout   ) ::  pts_ad    ! 1 : Adjoint of potential temperature [Celsius]
1833                                                                         ! 2 : Adjoint of salinity [psu]
1834      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   &
1835         & pn2_ad                   ! Adjoint of potential density (surface referenced)
1836      !
1837      !! * Local declarations
1838      INTEGER  ::   ji, jj, jk   ! dummy loop indices
1839      REAL(wp) ::             &
1840         zgde3w, zt, zs, zh,  &  ! temporary scalars
1841         zalbet, zbeta           !    "         "
1842      REAL(wp) ::             &
1843         ztad, zsad,          &  ! temporary scalars
1844         zalbetad, zbetaad       !    "         "
1845#if defined key_zdfddm
1846      REAL(wp) ::   zds, zdsad          ! temporary scalars
1847#endif
1848
1849      ! pn2_tl : interior points only (2=< jk =< jpkm1 )
1850      ! --------------------------
1851      zalbetad = 0.0_wp
1852      zbetaad  = 0.0_wp
1853      ztad     = 0.0_wp
1854      zsad     = 0.0_wp
1855#if defined key_zdfddm
1856      zdsad    = 0.0_wp
1857#endif
1858
1859      SELECT CASE ( nn_eos )
1860
1861      CASE ( 0 )               !== Jackett and McDougall (1994) formulation ==!
1862         DO jk = jpkm1, 2, -1
1863            DO jj = jpj, 1, -1
1864               DO ji = jpi, 1, -1
1865                  zgde3w = grav / fse3w(ji,jj,jk)
1866                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )          ! potential temperature at w-point
1867                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0   ! salinity anomaly (s-35) at w-point
1868                  zh = fsdepw(ji,jj,jk)                                     ! depth in meters  at w-point
1869
1870                  zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &   ! ratio alpha/beta
1871                     &                               - 0.203814e-03_wp ) * zt   &
1872                     &                               + 0.170907e-01_wp ) * zt   &
1873                     &   + 0.665157e-01_wp                                      &
1874                     &   +     ( - 0.678662e-05_wp * zs                         &
1875                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
1876                     &   +   ( ( - 0.302285e-13_wp * zh                         &
1877                     &           - 0.251520e-11_wp * zs                         &
1878                     &           + 0.512857e-12_wp * zt * zt           ) * zh   &
1879                     &           - 0.164759e-06_wp * zs                         &
1880                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
1881                     &                                  + 0.380374e-04_wp ) * zh
1882
1883                  zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta
1884                     &                               - 0.301985e-05_wp ) * zt      &
1885                     &   + 0.785567e-03_wp                                      &
1886                     &   + (     0.515032e-08_wp * zs                           &
1887                     &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     &
1888                     &   +(  (   0.121551e-17_wp * zh                           &
1889                     &         - 0.602281e-15_wp * zs                           &
1890                     &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     &
1891                     &                                + 0.408195e-10_wp   * zs     &
1892                     &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     &
1893                     &                                - 0.121555e-07_wp ) * zh
1894
1895#if defined key_zdfddm
1896
1897                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1898                  IF ( ABS( zds) <= 1.e-20 ) THEN
1899                     zds = 1.e-20
1900                     zdsad = 0.0_wp
1901                  ELSE
1902                     zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1903                     zdsad = rrau_ad(ji,jj,jk) * zalbet *( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds**2
1904                     zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1905                  ENDIF
1906                  pts_ad(ji,jj,jk-1,jp_tem) =  pts_ad(ji,jj,jk-1,jp_tem) +  rrau_ad(ji,jj,jk) * zalbet / zds
1907                  pts_ad(ji,jj,jk,jp_tem  ) =  pts_ad(ji,jj,jk,jp_tem  ) -  rrau_ad(ji,jj,jk) * zalbet / zds
1908                  zalbetad = zalbetad + rrau_ad(ji,jj,jk) * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1909                  rrau_ad(ji,jj,jk) = 0._wp
1910                  pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + zdsad
1911                  pts_ad(ji,jj,jk,jp_sal  ) = pts_ad(ji,jj,jk,jp_sal  ) - zdsad
1912                  zdsad = 0._wp
1913#endif
1914                     pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) + zalbet*zbeta*zgde3w*tmask(ji,jj,jk)*pn2_ad(ji,jj,jk)
1915                     pts_ad(ji,jj,jk,jp_tem  ) = pts_ad(ji,jj,jk,jp_tem  ) - zalbet*zbeta*zgde3w*tmask(ji,jj,jk)*pn2_ad(ji,jj,jk)
1916                     zalbetad = zalbetad + zbeta*zgde3w*tmask(ji,jj,jk)   &
1917                              &   *( pts  (ji,jj,jk-1,jp_tem) - pts  (ji,jj,jk,jp_tem) ) * pn2_ad(ji,jj,jk)
1918                     pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) - zbeta*tmask(ji,jj,jk)*zgde3w*pn2_ad(ji,jj,jk)
1919                     pts_ad(ji,jj,jk,jp_sal  ) = pts_ad(ji,jj,jk,jp_sal  ) + zbeta*tmask(ji,jj,jk)*zgde3w*pn2_ad(ji,jj,jk)
1920                     zbetaad = zbetaad &
1921                        & + zgde3w *tmask(ji,jj,jk)* (  zalbet * ( pts  (ji,jj,jk-1,jp_tem) - pts  (ji,jj,jk,jp_tem) )  &
1922                        & -  ( pts  (ji,jj,jk-1,jp_sal) - pts  (ji,jj,jk,jp_sal) ) )*pn2_ad(ji,jj,jk)
1923
1924                     pn2_ad(ji,jj,jk)  = 0.0_wp
1925
1926                     ztad = ztad + (    (     -3.*0.415613e-09_wp   * zt           &
1927                        &              +2.*0.555579e-07_wp ) * zt                  &
1928                        &      -           0.301985e-05_wp                         &
1929                        &      +           0.788212e-08_wp   * zs                  &
1930                        &      + (     -2.*0.213127e-11_wp   * zt                  &
1931                        &              -   0.175379e-14_wp   * zh                  &
1932                        &              +   0.192867e-09_wp         ) * zh ) *zbetaad
1933
1934                     zsad = zsad + (           2.*0.515032e-08_wp   * zs           &
1935                        &      +           0.788212e-08_wp   * zt                  &
1936                        &      -           0.356603e-06_wp                         &
1937                        &      + (     -   0.602281e-15_wp   * zh                  &
1938                        &              +   0.408195e-10_wp         ) * zh ) * zbetaad
1939
1940                     zbetaad = 0.0_wp
1941
1942                     ztad = ztad + (    (   ( -4.*0.255019e-07_wp   * zt           &! ratio alpha/beta
1943                        &              +3.*0.298357e-05_wp ) * zt                  &
1944                        &              -2.*0.203814e-03_wp ) * zt                  &
1945                        &      +           0.170907e-01_wp                         &
1946                        &      -           0.846960e-04_wp   * zs                  &
1947                        &      - (         0.933746e-06_wp                         &
1948                        &          - (  2.*0.791325e-08_wp                         &
1949                        &              +2.*0.512857e-12_wp   * zh ) * zt ) * zh    &
1950                        &                                                 ) *zalbetad
1951
1952                     zsad = zsad + (  -        2.*0.678662e-05_wp   * zs           &
1953                        &      -           0.846960e-04_wp   * zt                  &
1954                        &      +           0.378110e-02_wp                         &
1955                        &      + (     -   0.164759e-06_wp                         &
1956                        &              -   0.251520e-11_wp   * zh ) * zh           &
1957                        &                                                 ) *zalbetad
1958
1959                     zalbetad = 0.0_wp
1960
1961
1962                     pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + 0.5 * zsad
1963                     pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + 0.5 * zsad
1964                     zsad = 0.0_wp
1965
1966                     pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + 0.5 * ztad
1967                     pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) + 0.5 * ztad
1968                     ztad = 0.0_wp
1969
1970               END DO
1971            END DO
1972         END DO
1973         !
1974      CASE ( 1 )               !==  Linear formulation = F( temperature )  ==!
1975         DO jk = jpkm1, 2, -1
1976            pts_ad(:,:,jk-1,jp_tem) = pts_ad(:,:,jk-1,jp_tem)  + rn_alpha * pn2_ad(:,:,jk) &
1977                                      & * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1978            pts_ad(:,:,jk,jp_tem  ) = pts_ad(:,:,jk,jp_tem  )  - rn_alpha * pn2_ad(:,:,jk) &
1979                                      & * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1980            pn2_ad(:,:,jk) = 0.0_wp
1981         END DO
1982         !
1983      CASE ( 2 )               !==  Linear formulation = F( temperature , salinity )  ==!
1984#if defined key_zdfddm
1985         DO jk = jpkm1, 2, -1
1986            DO jj = jpj, 1, -1
1987               DO ji = jpi, 1, -1
1988                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1989                  IF ( ABS( zds) <= 1.e-20 ) THEN
1990                     zds = 1.e-20
1991                     zdsad = 0.0_wp
1992                  ELSE
1993                     zdsad = zdsad - rrau_ad(ji,jj,jk) * ralpbet &
1994                     & * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds**2
1995                  ENDIF
1996                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1997                  pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem)   &
1998                     & + rrau_ad(ji,jj,jk) * ralpbet / zds
1999                  pts_ad(ji,jj,jk,jp_tem  ) = pts_ad(ji,jj,jk,jp_tem  )   &
2000                     & - rrau_ad(ji,jj,jk) * ralpbet / zds
2001                  rrau_ad(ji,jj,jk)  = 0._wp
2002
2003                  pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + zdsad
2004                  pts_ad(ji,jj,jk,jp_sal  ) = pts_ad(ji,jj,jk,jp_sal  ) - zdsad
2005                  zdsad = 0._wp
2006               END DO
2007            END DO
2008         END DO
2009#endif
2010         DO jk = jpkm1, 2, -1
2011            pts_ad(:,:,jk-1,jp_tem) = pts_ad(:,:,jk-1,jp_tem) + rn_alpha * pn2_ad(:,:,jk) &
2012                       & * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2013            pts_ad(:,:,jk,jp_tem  ) = pts_ad(:,:,jk,jp_tem  ) - rn_alpha * pn2_ad(:,:,jk) &
2014                       & * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2015            pts_ad(:,:,jk-1,jp_sal) = pts_ad(:,:,jk-1,jp_sal) - rn_beta  * pn2_ad(:,:,jk) &
2016                       & * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2017            pts_ad(:,:,jk,jp_sal  ) = pts_ad(:,:,jk,jp_sal  ) + rn_beta  * pn2_ad(:,:,jk) &
2018                       & * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2019            pn2_ad(:,:,jk) = 0.0_wp
2020         END DO
2021      END SELECT
2022   END SUBROUTINE eos_bn2_adj
2023
2024   SUBROUTINE eos_alpbet_tan( pts, pts_tl, palpbet_tl, beta0_tl )
2025      !!----------------------------------------------------------------------
2026      !!                 ***  ROUTINE eos_alpbet_tan  ***
2027      !!
2028      !! ** Purpose of the direct routine :
2029      !!       Calculates the in situ thermal/haline expansion ratio at T-points
2030      !!
2031      !! ** Method  of the direct routine :
2032      !!       calculates alpha / beta ratio at T-points
2033      !!       * nn_eos = 0  : UNESCO sea water properties
2034      !!                       The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
2035      !!                       polynomial expression of McDougall (1987).
2036      !!                       Scalar beta0 is returned = 1.
2037      !!       * nn_eos = 1  : linear equation of state (temperature only)
2038      !!                       The ratio is undefined, so we return alpha as palpbet
2039      !!                       Scalar beta0 is returned = 0.
2040      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
2041      !!                       The alpha/beta ratio is returned as ralpbet
2042      !!                       Scalar beta0 is returned = 1.
2043      !!
2044      !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points
2045      !!            :   beta0   : 1. or 0.
2046      !!----------------------------------------------------------------------
2047      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts_tl, &    ! linear tangent of pot. temperature & salinity
2048         &                                                      pts          ! pot. temperature & salinity
2049      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpbet_tl   ! thermal/haline expansion ratio
2050      REAL(wp),                              INTENT(  out) ::   beta0_tl     ! set = 1 except with case 1 eos, rho=rho(T)
2051      !!
2052      INTEGER  ::   ji, jj, jk      ! dummy loop indices
2053      REAL(wp) ::   zt, zs, zh, &   ! local scalars
2054         &          zttl, zstl
2055      !!----------------------------------------------------------------------
2056      !
2057      IF( nn_timing == 1 ) CALL timing_start('eos_alpbet_tan')
2058      !
2059      SELECT CASE ( nn_eos )
2060      !
2061      CASE ( 0 )               ! Jackett and McDougall (1994) formulation
2062         DO jk = 1, jpk
2063            DO jj = 1, jpj
2064               DO ji = 1, jpi
2065                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature
2066                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35)
2067                  zh = fsdept(ji,jj,jk)               ! depth in meters
2068                  !! Tangent part
2069                  zttl = pts_tl(ji,jj,jk,jp_tem)           ! potential temperature
2070                  zstl = pts_tl(ji,jj,jk,jp_sal)           ! salinity anomaly (s-35)
2071                  palpbet_tl(ji,jj,jk) =                                           &
2072                     &     ( ( ( ( - 4. * 0.255019e-07_wp * zt                     &
2073                     &           + 3. * 0.298357e-05_wp ) * zt                     &
2074                     &           - 2. * 0.203814e-03_wp ) * zt                     &
2075                     &           + 0.170907e-01_wp * zt )                          &
2076                     &           - 0.846960e-04_wp * zs                            &
2077                     &     + ( ( 2. * 0.512857e-12_wp * zt ) * zh                  &
2078                     &     +   ( 2. * 0.791325e-08_wp * zt                         &
2079                     &           - 0.933746e-06_wp ) ) * zh ) * zttl               &
2080                     &     + ( ( - 2. * 0.678662e-05_wp * zs                       &
2081                     &           - 0.846960e-04_wp * zt                            &
2082                     &           + 0.378110e-02_wp )                               &
2083                     &       + ( - 0.251520e-11_wp  * zh                           &
2084                     &           - 0.164759e-06_wp ) * zh ) * zstl
2085               END DO
2086            END DO
2087         END DO
2088         beta0_tl = 0._wp
2089         !
2090      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==!
2091         palpbet_tl(:,:,:) = 0._wp
2092         beta0_tl = 0._wp
2093         !
2094      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==!
2095         palpbet_tl(:,:,:) = 0._wp
2096         beta0_tl = 0._wp
2097         !
2098      CASE DEFAULT
2099         IF(lwp) WRITE(numout,cform_err)
2100         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
2101         nstop = nstop + 1
2102         !
2103      END SELECT
2104      !
2105      IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet_tan')
2106      !
2107   END SUBROUTINE eos_alpbet_tan
2108
2109   SUBROUTINE eos_alpbet_adj( pts, pts_ad, palpbet_ad, beta0_ad )
2110      !!----------------------------------------------------------------------
2111      !!                 ***  ROUTINE eos_alpbet_adj  ***
2112      !!
2113      !! ** Purpose of the direct routine :
2114      !!       Calculates the in situ thermal/haline expansion ratio at T-points
2115      !!
2116      !! ** Method  of the direct routine :
2117      !!       calculates alpha / beta ratio at T-points
2118      !!       * nn_eos = 0  : UNESCO sea water properties
2119      !!                       The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
2120      !!                       polynomial expression of McDougall (1987).
2121      !!                       Scalar beta0 is returned = 1.
2122      !!       * nn_eos = 1  : linear equation of state (temperature only)
2123      !!                       The ratio is undefined, so we return alpha as palpbet
2124      !!                       Scalar beta0 is returned = 0.
2125      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
2126      !!                       The alpha/beta ratio is returned as ralpbet
2127      !!                       Scalar beta0 is returned = 1.
2128      !!
2129      !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points
2130      !!            :   beta0   : 1. or 0.
2131      !!----------------------------------------------------------------------
2132      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) ::   pts_ad       ! linear tangent of pot. temperature & salinity
2133      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in)    ::   pts          ! pot. temperature & salinity
2134      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(inout) ::   palpbet_ad   ! thermal/haline expansion ratio
2135      REAL(wp),                              INTENT(inout) ::   beta0_ad     ! set = 1 except with case 1 eos, rho=rho(T)
2136      !!
2137      INTEGER  ::   ji, jj, jk      ! dummy loop indices
2138      REAL(wp) ::   zt, zs, zh, &   ! local scalars
2139         &          ztad, zsad
2140      !!----------------------------------------------------------------------
2141      !
2142      ztad = 0.0_wp
2143      zsad = 0.0_wp
2144      IF( nn_timing == 1 ) CALL timing_start('eos_alpbet_adj')
2145      !
2146      SELECT CASE ( nn_eos )
2147      !
2148      CASE ( 0 )               ! Jackett and McDougall (1994) formulation
2149         DO jk = jpk, 1, -1
2150            DO jj = jpj, 1, -1
2151               DO ji = jpi, 1, -1
2152                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature
2153                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35)
2154                  zh = fsdept(ji,jj,jk)               ! depth in meters
2155                  !! Adjoint part
2156                  ztad = ztad + ( ( ( ( - 4. * 0.255019e-07_wp * zt                &
2157                     &           + 3. * 0.298357e-05_wp ) * zt                     &
2158                     &           - 2. * 0.203814e-03_wp ) * zt                     &
2159                     &           + 0.170907e-01_wp * zt )                          &
2160                     &           - 0.846960e-04_wp * zs                            &
2161                     &     + ( ( 2. * 0.512857e-12_wp * zt ) * zh                  &
2162                     &     +   ( 2. * 0.791325e-08_wp * zt                         &
2163                     &           - 0.933746e-06_wp ) ) * zh ) * palpbet_ad(ji,jj,jk)
2164                  zsad = zsad + ( ( - 2. * 0.678662e-05_wp * zs                    &
2165                     &           - 0.846960e-04_wp * zt                            &
2166                     &           + 0.378110e-02_wp )                               &
2167                     &       + ( - 0.251520e-11_wp  * zh                           &
2168                     &           - 0.164759e-06_wp ) * zh ) * palpbet_ad(ji,jj,jk)
2169                  palpbet_ad(ji,jj,jk) = 0.0_wp
2170                  pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + ztad
2171                  pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + zsad
2172                  ztad = 0.0_wp
2173                  zsad = 0.0_wp
2174               END DO
2175            END DO
2176         END DO
2177         beta0_ad = 0._wp
2178         !
2179      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==!
2180         palpbet_ad(:,:,:) = 0._wp
2181         beta0_ad = 0._wp
2182         !
2183      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==!
2184         palpbet_ad(:,:,:) = 0._wp
2185         beta0_ad = 0._wp
2186         !
2187      CASE DEFAULT
2188         IF(lwp) WRITE(numout,cform_err)
2189         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
2190         nstop = nstop + 1
2191         !
2192      END SELECT
2193      !
2194      IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet_adj')
2195      !
2196   END SUBROUTINE eos_alpbet_adj
2197
2198#if defined key_tam
2199   SUBROUTINE eos_insitu_adj_tst( kumadt )
2200      !!-----------------------------------------------------------------------
2201      !!
2202      !!                  ***  ROUTINE eos_adj_tst ***
2203      !!
2204      !! ** Purpose : Test the adjoint routine.
2205      !!
2206      !! ** Method  : Verify the scalar product
2207      !!
2208      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
2209      !!
2210      !!              where  L   = tangent routine
2211      !!                     L^T = adjoint routine
2212      !!                     W   = diagonal matrix of scale factors
2213      !!                     dx  = input perturbation (random field)
2214      !!                     dy  = L dx
2215      !!
2216      !!
2217      !! History :
2218      !!        ! 08-07 (A. Vidard)
2219      !!-----------------------------------------------------------------------
2220      !! * Modules used
2221
2222      !! * Arguments
2223      INTEGER, INTENT(IN) :: &
2224         & kumadt             ! Output unit
2225      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2226         zts                       ! potential temperature
2227                                   ! salinity
2228      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2229         & zts_adout                 ! potential temperature
2230                                    ! salinity
2231      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2232         & zrd_adin                 ! anomaly density
2233      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2234         & zts_tlin                 ! potential temperature
2235                                    ! salinity
2236      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2237         & znts                     ! potential temperature
2238                                   ! salinity
2239      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2240         & zrd_tlout                  ! anomaly density
2241      REAL(KIND=wp) ::       &
2242         & zsp1,             & ! scalar product involving the tangent routine
2243         & zsp2                ! scalar product involving the adjoint routine
2244      INTEGER :: &
2245         & ji, &
2246         & jj, &
2247         & jk, &
2248      & jn, &
2249      & jeos
2250      CHARACTER(LEN=14) :: cl_name
2251      ALLOCATE( &
2252         & zts(     jpi, jpj, jpk, jpts ),  &
2253         & znts(      jpi, jpj, jpk, jpts ),  &
2254         & zts_adout( jpi, jpj, jpk,jpts ),  &
2255         & zrd_adin( jpi, jpj, jpk ),  &
2256         & zts_tlin(  jpi, jpj, jpk,jpts ),  &
2257         & zrd_tlout(jpi, jpj, jpk )    )
2258      ! Initialize the reference state
2259      zts = tsn
2260      ! store initial nn_eos
2261      jeos = nn_eos
2262      DO jn = 0, 2
2263         nn_eos = jn
2264         !=============================================================
2265         ! 1) dx = ( T ) and dy = ( T )
2266         !=============================================================
2267
2268         !--------------------------------------------------------------------
2269         ! Reset the tangent and adjoint variables
2270         !--------------------------------------------------------------------
2271         zts_tlin(:,:,:,:)     = 0.0_wp
2272         zrd_tlout(:,:,:)   = 0.0_wp
2273         zts_adout(:,:,:,:)    = 0.0_wp
2274         zrd_adin(:,:,:)    = 0.0_wp
2275
2276         !--------------------------------------------------------------------
2277         ! Initialize the tangent input with random noise: dx
2278         !--------------------------------------------------------------------
2279         CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2280         CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2281
2282         DO jk = 1, jpk
2283            DO jj = nldj, nlej
2284               DO ji = nldi, nlei
2285                  zts_tlin(ji,jj,jk,jp_tem) = znts(ji,jj,jk,jp_tem)
2286                  zts_tlin(ji,jj,jk,jp_sal) = znts(ji,jj,jk,jp_sal)
2287               END DO
2288            END DO
2289         END DO
2290         CALL eos_insitu_tan(zts, zts_tlin, zrd_tlout)
2291
2292         DO jk = 1, jpk
2293            DO jj = nldj, nlej
2294               DO ji = nldi, nlei
2295                  zrd_adin(ji,jj,jk)   = zrd_tlout(ji,jj,jk) &
2296                       &                 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2297                       &                 * tmask(ji,jj,jk)
2298               END DO
2299            END DO
2300         END DO
2301
2302         !--------------------------------------------------------------------
2303         ! Compute the scalar product: ( L dx )^T W dy
2304         !--------------------------------------------------------------------
2305
2306         zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2307
2308         !--------------------------------------------------------------------
2309         ! Call the adjoint routine: dx^* = L^T dy^*
2310         !--------------------------------------------------------------------
2311         CALL eos_insitu_adj(zts, zts_adout, zrd_adin)
2312         zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2313                & DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2314
2315         ! Compare the scalar products
2316
2317         ! Compare the scalar products
2318         ! 14 char:'12345678901234'
2319         SELECT CASE( jn )
2320         CASE (0) ; cl_name = 'eos_adj ins T1'
2321         CASE (1) ; cl_name = 'eos_adj ins T2'
2322         CASE (2) ; cl_name = 'eos_adj ins T3'
2323         END SELECT
2324         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2325      ENDDO
2326      ! restore initial nn_eos
2327      nn_eos = jeos
2328
2329      ! Deallocate memory
2330
2331      DEALLOCATE(      &
2332         & zts,       &
2333         & zts_adout,   &
2334         & zrd_adin,   &
2335         & zts_tlin,    &
2336         & zrd_tlout,  &
2337         & znts        &
2338         & )
2339   END SUBROUTINE eos_insitu_adj_tst
2340
2341   SUBROUTINE eos_insitu_pot_adj_tst( kumadt )
2342      !!-----------------------------------------------------------------------
2343      !!
2344      !!                  ***  ROUTINE eos_adj_tst ***
2345      !!
2346      !! ** Purpose : Test the adjoint routine.
2347      !!
2348      !! ** Method  : Verify the scalar product
2349      !!
2350      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
2351      !!
2352      !!              where  L   = tangent routine
2353      !!                     L^T = adjoint routine
2354      !!                     W   = diagonal matrix of scale factors
2355      !!                     dx  = input perturbation (random field)
2356      !!                     dy  = L dx
2357      !!
2358      !! ** Action  : Separate tests are applied for the following dx and dy:
2359      !!
2360      !!              1) dx = ( SSH ) and dy = ( SSH )
2361      !!
2362      !! History :
2363      !!        ! 08-07 (A. Vidard)
2364      !!-----------------------------------------------------------------------
2365      !! * Modules used
2366
2367      !! * Arguments
2368      INTEGER, INTENT(IN) :: &
2369         & kumadt             ! Output unit
2370      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2371         zts                        ! potential temperature
2372                                    ! salinity
2373      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2374         & zts_adout                ! potential temperature
2375                                    ! salinity
2376      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2377         & zrd_adin                 ! anomaly density
2378      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2379         & zrhop_adin               ! volume mass
2380      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2381         & zts_tlin                  ! potential temperature
2382                                     ! salinity
2383      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2384         & znts                     ! potential temperature
2385                                    ! salinity
2386      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2387         & zrd_tlout                  ! anomaly density
2388      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2389         & zrhop_tlout                  ! volume mass
2390      REAL(KIND=wp) ::       &
2391         & zsp1,             & ! scalar product involving the tangent routine
2392         & zsp2                ! scalar product involving the adjoint routine
2393      INTEGER :: &
2394         & ji, &
2395         & jj, &
2396         & jk, &
2397         & jn, &
2398         & jeos
2399      CHARACTER(LEN=14) :: cl_name
2400
2401       ! Allocate memory
2402      ALLOCATE( &
2403         & zts(     jpi, jpj, jpk, jpts ),  &
2404         & zts_adout( jpi, jpj, jpk, jpts ),  &
2405         & zrhop_adin( jpi, jpj, jpk ),  &
2406         & zrd_adin( jpi, jpj, jpk ),  &
2407         & zts_tlin(  jpi, jpj, jpk, jpts ),  &
2408         & znts(      jpi, jpj, jpk, jpts ),  &
2409         & zrd_tlout(jpi, jpj, jpk ),  &
2410         & zrhop_tlout(jpi, jpj, jpk )    )
2411
2412      ! Initialize random field standard deviationsthe reference state
2413      zts = tsn
2414
2415      ! store initial nn_eos
2416      jeos = nn_eos
2417      DO jn = 0, 2
2418         nn_eos = jn
2419         !=============================================
2420         !  testing of eos_insitu_pot
2421         !=============================================
2422
2423         !=============================================================
2424         ! 1) dx = ( T ) and dy = ( T )
2425         !=============================================================
2426
2427         !--------------------------------------------------------------------
2428         ! Reset the tangent and adjoint variables
2429         !--------------------------------------------------------------------
2430         zts_tlin(:,:,:,:)     = 0.0_wp
2431         zrd_tlout(:,:,:)   = 0.0_wp
2432         zrhop_tlout(:,:,:) = 0.0_wp
2433         zts_adout(:,:,:,:)    = 0.0_wp
2434         zrhop_adin(:,:,:)  = 0.0_wp
2435         zrd_adin(:,:,:)    = 0.0_wp
2436
2437         !--------------------------------------------------------------------
2438         ! Initialize the tangent input with random noise: dx
2439         !--------------------------------------------------------------------
2440         CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2441         CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2442         DO jk = 1, jpk
2443            DO jj = nldj, nlej
2444               DO ji = nldi, nlei
2445                  zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2446               END DO
2447            END DO
2448         END DO
2449         !--------------------------------------------------------------------
2450         ! Call the tangent routine: dy = L dx
2451         !--------------------------------------------------------------------
2452
2453         call eos_insitu_pot_tan ( zts, zts_tlin, zrd_tlout, zrhop_tlout )
2454
2455         !--------------------------------------------------------------------
2456         ! Initialize the adjoint variables: dy^* = W dy
2457         !--------------------------------------------------------------------
2458         DO jk = 1, jpk
2459            DO jj = nldj, nlej
2460               DO ji = nldi, nlei
2461                  zrd_adin(ji,jj,jk)   = zrd_tlout(ji,jj,jk) &
2462                       &                 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2463                       &                 * tmask(ji,jj,jk)
2464                  zrhop_adin(ji,jj,jk) = zrhop_tlout(ji,jj,jk) &
2465                       &                 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2466                       &                 * tmask(ji,jj,jk)
2467               END DO
2468            END DO
2469         END DO
2470
2471         !--------------------------------------------------------------------
2472         ! Compute the scalar product: ( L dx )^T W dy
2473         !--------------------------------------------------------------------
2474
2475         zsp1 =  DOT_PRODUCT( zrd_tlout  , zrd_adin   ) &
2476              &  + DOT_PRODUCT( zrhop_tlout, zrhop_adin )
2477         !--------------------------------------------------------------------
2478         ! Call the adjoint routine: dx^* = L^T dy^*
2479         !--------------------------------------------------------------------
2480
2481         CALL eos_insitu_pot_adj( zts, zts_adout, zrd_adin, zrhop_adin )
2482         !--------------------------------------------------------------------
2483         ! Compute the scalar product: dx^T L^T W dy
2484         !--------------------------------------------------------------------
2485
2486         zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2487                & DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2488         ! Compare the scalar products
2489
2490         ! 14 char:'12345678901234'
2491         SELECT CASE( jn )
2492         CASE (0) ; cl_name = 'eos_adj pot T1'
2493         CASE (1) ; cl_name = 'eos_adj pot T2'
2494         CASE (2) ; cl_name = 'eos_adj pot T3'
2495         END SELECT
2496         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2497
2498      ENDDO
2499
2500      ! restore initial nn_eos
2501      nn_eos = jeos
2502
2503      ! Deallocate memory
2504      DEALLOCATE(       &
2505         & zts,       &
2506         & zts_adout,   &
2507         & zrd_adin,   &
2508         & zrhop_adin, &
2509         & zts_tlin,    &
2510         & zrd_tlout,  &
2511         & zrhop_tlout,&
2512         & znts    )
2513   END SUBROUTINE eos_insitu_pot_adj_tst
2514
2515   SUBROUTINE eos_alpbet_adj_tst(kumadt)
2516      !!-----------------------------------------------------------------------
2517      !!
2518      !!                  ***  ROUTINE eos_adj_tst ***
2519      !!
2520      !! ** Purpose : Test the adjoint routine.
2521      !!
2522      !! ** Method  : Verify the scalar product
2523      !!
2524      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
2525      !!
2526      !!              where  L   = tangent routine
2527      !!                     L^T = adjoint routine
2528      !!                     W   = diagonal matrix of scale factors
2529      !!                     dx  = input perturbation (random field)
2530      !!                     dy  = L dx
2531      !!
2532      !! ** Action  : Separate tests are applied for the following dx and dy:
2533      !!
2534      !!              1) dx = ( SSH ) and dy = ( SSH )
2535      !!
2536      !! History :
2537      !!        ! 05-12 (P.-A. Bouttier)
2538      !!-----------------------------------------------------------------------
2539      !! * Modules used
2540
2541      !! * Arguments
2542      INTEGER, INTENT(IN) :: &
2543         & kumadt             ! Output unit
2544      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2545         & zts                        ! potential temperature
2546                                    ! salinity
2547      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2548         & zts_adout                ! potential temperature
2549                                    ! salinity
2550      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2551         & zpalpbet_adin
2552      REAL(wp)                       ::   &
2553         & zbeta0_adin
2554      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2555         & zts_tlin                  ! potential temperature
2556                                     ! salinity
2557      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2558         & znts                     ! potential temperature
2559                                    ! salinity
2560      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2561         & zpalpbet_tlout
2562      REAL(wp)                       ::   &
2563         & zbeta0_tlout
2564      REAL(KIND=wp) ::       &
2565         & zsp1,             & ! scalar product involving the tangent routine
2566         & zsp2                ! scalar product involving the adjoint routine
2567      INTEGER :: &
2568         & ji, &
2569         & jj, &
2570         & jk, &
2571         & jn, &
2572         & jeos
2573      CHARACTER(LEN=14) :: cl_name
2574
2575       ! Allocate memory
2576      ALLOCATE( &
2577         & zts(     jpi, jpj, jpk, jpts ),  &
2578         & zts_adout( jpi, jpj, jpk, jpts ),  &
2579         & zpalpbet_adin( jpi, jpj, jpk ),  &
2580         & zts_tlin(  jpi, jpj, jpk, jpts ),  &
2581         & znts(      jpi, jpj, jpk, jpts ),  &
2582         & zpalpbet_tlout(jpi, jpj, jpk ) )
2583
2584
2585      ! Initialize random field standard deviationsthe reference state
2586      zts = tsn
2587
2588      ! store initial nn_eos
2589      jeos = nn_eos
2590      DO jn = 0, 2
2591         nn_eos = jn
2592         !=============================================
2593         !  testing of eos_insitu_pot
2594         !=============================================
2595
2596         !=============================================================
2597         ! 1) dx = ( T ) and dy = ( T )
2598         !=============================================================
2599
2600         !--------------------------------------------------------------------
2601         ! Reset the tangent and adjoint variables
2602         !--------------------------------------------------------------------
2603         zts_tlin(:,:,:,:)     = 0.0_wp
2604         zpalpbet_tlout(:,:,:)   = 0.0_wp
2605         zbeta0_tlout = 0.0_wp
2606         zts_adout(:,:,:,:)    = 0.0_wp
2607         zbeta0_adin  = 0.0_wp
2608         zpalpbet_adin(:,:,:)    = 0.0_wp
2609
2610         !--------------------------------------------------------------------
2611         ! Initialize the tangent input with random noise: dx
2612         !--------------------------------------------------------------------
2613         CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2614         CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2615         DO jk = 1, jpk
2616            DO jj = nldj, nlej
2617               DO ji = nldi, nlei
2618                  zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2619               END DO
2620            END DO
2621         END DO
2622         !--------------------------------------------------------------------
2623         ! Call the tangent routine: dy = L dx
2624         !--------------------------------------------------------------------
2625
2626         call eos_alpbet_tan ( zts, zts_tlin, zpalpbet_tlout, zbeta0_tlout )
2627
2628         !--------------------------------------------------------------------
2629         ! Initialize the adjoint variables: dy^* = W dy
2630         !--------------------------------------------------------------------
2631         DO jk = 1, jpk
2632            DO jj = nldj, nlej
2633               DO ji = nldi, nlei
2634                  zpalpbet_adin(ji,jj,jk)   = zpalpbet_tlout(ji,jj,jk) &
2635                       &                 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2636                       &                 * tmask(ji,jj,jk)
2637               END DO
2638            END DO
2639         END DO
2640         zbeta0_adin = zbeta0_tlout
2641
2642         !--------------------------------------------------------------------
2643         ! Compute the scalar product: ( L dx )^T W dy
2644         !--------------------------------------------------------------------
2645
2646         zsp1 =  DOT_PRODUCT( zpalpbet_tlout  , zpalpbet_adin   ) &
2647              &  + zbeta0_tlout *  zbeta0_adin
2648         !--------------------------------------------------------------------
2649         ! Call the adjoint routine: dx^* = L^T dy^*
2650         !--------------------------------------------------------------------
2651
2652         CALL eos_alpbet_adj( zts, zts_adout, zpalpbet_adin, zbeta0_adin )
2653
2654         !--------------------------------------------------------------------
2655         ! Compute the scalar product: dx^T L^T W dy
2656         !--------------------------------------------------------------------
2657
2658         zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2659                & DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2660         ! Compare the scalar products
2661
2662         ! 14 char:'12345678901234'
2663         SELECT CASE( jn )
2664         CASE (0) ; cl_name = 'eos_adj ab  T1'
2665         CASE (1) ; cl_name = 'eos_adj ab  T2'
2666         CASE (2) ; cl_name = 'eos_adj ab  T3'
2667         END SELECT
2668         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2669
2670      ENDDO
2671
2672      ! restore initial nn_eos
2673      nn_eos = jeos
2674
2675      ! Deallocate memory
2676      DEALLOCATE(       &
2677         & zts,       &
2678         & zts_adout,   &
2679         & zpalpbet_adin,   &
2680         & zts_tlin,    &
2681         & zpalpbet_tlout,  &
2682         & znts    )
2683
2684   END SUBROUTINE eos_alpbet_adj_tst
2685
2686   SUBROUTINE eos_insitu_2d_adj_tst( kumadt )
2687      !!-----------------------------------------------------------------------
2688      !!
2689      !!                  ***  ROUTINE eos_adj_tst ***
2690      !!
2691      !! ** Purpose : Test the adjoint routine.
2692      !!
2693      !! ** Method  : Verify the scalar product
2694      !!
2695      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
2696      !!
2697      !!              where  L   = tangent routine
2698      !!                     L^T = adjoint routine
2699      !!                     W   = diagonal matrix of scale factors
2700      !!                     dx  = input perturbation (random field)
2701      !!                     dy  = L dx
2702      !!
2703      !! ** Action  : Separate tests are applied for the following dx and dy:
2704      !!
2705      !!              1) dx = ( SSH ) and dy = ( SSH )
2706      !!
2707      !! History :
2708      !!        ! 08-07 (A. Vidard)
2709      !!-----------------------------------------------------------------------
2710      !! * Modules used
2711
2712      !! * Arguments
2713      INTEGER, INTENT(IN) :: &
2714         & kumadt             ! Output unit
2715      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   &
2716         zdep                       ! depth
2717      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2718         zts                        ! potential temperature
2719                                    ! salinity
2720      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2721         & zts_adout                ! potential temperature
2722                                    ! salinity
2723      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   &
2724         & zrd_adin                 ! anomaly density
2725      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2726         & zts_tlin                 ! potential temperature
2727                                    ! salinity
2728      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2729         & znts                     ! potential temperature
2730                                    ! salinity
2731      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   &
2732         & zrd_tlout                  ! anomaly density
2733      REAL(KIND=wp) ::       &
2734         & zsp1,             & ! scalar product involving the tangent routine
2735         & zsp2                ! scalar product involving the adjoint routine
2736      INTEGER :: &
2737         & ji, &
2738         & jj, &
2739         & jn, &
2740         & jeos
2741      CHARACTER(LEN=14) :: cl_name
2742      ! Allocate memory
2743
2744      ALLOCATE( &
2745         & zdep(     jpi, jpj),  &
2746         & zts(     jpi, jpj, jpts ),  &
2747         & znts(      jpi, jpj, jpts ),  &
2748         & zts_adout( jpi, jpj, jpts ),  &
2749         & zrd_adin( jpi, jpj ),  &
2750         & zts_tlin(  jpi, jpj, jpts ),  &
2751         & zrd_tlout(jpi, jpj )    )
2752
2753      ! Initialize the reference state
2754      zts(:,:,:) = tsn(:,:,2,:)
2755      zdep(:,:) = fsdept(:,:,2)
2756
2757      ! store initial nn_eos
2758      jeos = nn_eos
2759      DO jn = 0, 2
2760         nn_eos = jn
2761         !=============================================================
2762         ! 1) dx = ( T ) and dy = ( T )
2763         !=============================================================
2764
2765         !--------------------------------------------------------------------
2766         ! Reset the tangent and adjoint variables
2767         !--------------------------------------------------------------------
2768         zts_tlin(:,:,:)     = 0.0_wp
2769         zrd_tlout(:,:)   = 0.0_wp
2770         zts_adout(:,:,:)    = 0.0_wp
2771         zrd_adin(:,:)    = 0.0_wp
2772
2773         !--------------------------------------------------------------------
2774         ! Initialize the tangent input with random noise: dx
2775         !--------------------------------------------------------------------
2776         CALL grid_random( znts(:,:,jp_tem), 'T', 0.0_wp, stdt )
2777         CALL grid_random( znts(:,:,jp_sal), 'T', 0.0_wp, stds )
2778         DO jj = nldj, nlej
2779            DO ji = nldi, nlei
2780               zts_tlin(ji,jj,:) = znts(ji,jj,:)
2781            END DO
2782         END DO
2783
2784         CALL eos_insitu_2d_tan(zts, zdep, zts_tlin, zrd_tlout)
2785
2786         DO jj = nldj, nlej
2787            DO ji = nldi, nlei
2788               zrd_adin(ji,jj)   = zrd_tlout(ji,jj) &
2789                    &                 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,2)&
2790                    &                 * tmask(ji,jj,2)
2791            END DO
2792         END DO
2793
2794         !--------------------------------------------------------------------
2795         ! Compute the scalar product: ( L dx )^T W dy
2796         !--------------------------------------------------------------------
2797
2798         zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2799
2800         !--------------------------------------------------------------------
2801         ! Call the adjoint routine: dx^* = L^T dy^*
2802         !--------------------------------------------------------------------
2803         CALL eos_insitu_2d_adj(zts, zdep, zts_adout, zrd_adin)
2804         zsp2 = DOT_PRODUCT( zts_tlin(:,:,jp_tem), zts_adout(:,:,jp_tem) ) + &
2805              & DOT_PRODUCT( zts_tlin(:,:,jp_sal), zts_adout(:,:,jp_sal) )
2806
2807         ! Compare the scalar products
2808
2809         ! 14 char:'12345678901234'
2810         SELECT CASE( jn )
2811         CASE (0) ; cl_name = 'eos_adj 2d  T1'
2812         CASE (1) ; cl_name = 'eos_adj 2d  T2'
2813         CASE (2) ; cl_name = 'eos_adj 2d  T3'
2814         END SELECT
2815         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2816
2817      ENDDO
2818
2819      ! restore initial nn_eos
2820      nn_eos = jeos
2821
2822      ! Deallocate memory
2823
2824      DEALLOCATE(      &
2825         & zdep,       &
2826         & zts,       &
2827         & zts_adout,   &
2828         & zrd_adin,   &
2829         & zts_tlin,    &
2830         & zrd_tlout,  &
2831         & znts   )
2832   END SUBROUTINE eos_insitu_2d_adj_tst
2833
2834   SUBROUTINE eos_adj_tst( kumadt )
2835      !!-----------------------------------------------------------------------
2836      !!
2837      !!                  ***  ROUTINE eos_adj_tst ***
2838      !!
2839      !! ** Purpose : Test the adjoint routine.
2840      !!
2841      !! History :
2842      !!        ! 08-07 (A. Vidard)
2843      !!-----------------------------------------------------------------------
2844      !! * Arguments
2845      INTEGER, INTENT(IN) :: &
2846         & kumadt             ! Output unit
2847
2848      CALL eos_insitu_adj_tst( kumadt )
2849      CALL eos_insitu_pot_adj_tst( kumadt )
2850      CALL eos_insitu_2d_adj_tst( kumadt )
2851      CALL eos_alpbet_adj_tst( kumadt )
2852
2853   END SUBROUTINE eos_adj_tst
2854
2855   SUBROUTINE bn2_adj_tst( kumadt )
2856      !!-----------------------------------------------------------------------
2857      !!
2858      !!                  ***  ROUTINE bn2_adj_tst ***
2859      !!
2860      !! ** Purpose : Test the adjoint routine.
2861      !!
2862      !! ** Method  : Verify the scalar product
2863      !!
2864      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
2865      !!
2866      !!              where  L   = tangent routine
2867      !!                     L^T = adjoint routine
2868      !!                     W   = diagonal matrix of scale factors
2869      !!                     dx  = input perturbation (random field)
2870      !!                     dy  = L dx
2871      !!
2872      !! ** Action  : Separate tests are applied for the following dx and dy:
2873      !!
2874      !!              1) dx = ( SSH ) and dy = ( SSH )
2875      !!
2876      !! History :
2877      !!        ! 08-07 (A. Vidard)
2878      !!-----------------------------------------------------------------------
2879      !! * Modules used
2880
2881      !! * Arguments
2882      INTEGER, INTENT(IN) :: &
2883         & kumadt             ! Output unit
2884      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2885         zts                      ! potential temperature
2886                                  ! salinity
2887      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2888         & zts_adout              ! potential temperature
2889                                  ! salinity
2890      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2891         & zrd_adin,               &  ! potential density (surface referenced)
2892         & zrd_adout                  ! potential density (surface referenced)
2893      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2894         & zts_tlin,              &  ! potential temperature
2895                                     ! salinity
2896         & zts_tlout                 ! potential temperature
2897                                     ! salinity
2898      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   &
2899         & zrd_tlout                  ! potential density (surface referenced)
2900      REAL(KIND=wp) ::       &
2901         & zsp1,             & ! scalar product involving the tangent routine
2902         & zsp2                ! scalar product involving the adjoint routine
2903      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   &
2904         & znts                   ! potential temperature
2905                                  ! salinity
2906      INTEGER :: &
2907         & ji, &
2908         & jj, &
2909         & jk, &
2910         & jn, &
2911         & jeos
2912      CHARACTER(LEN=14) :: cl_name
2913
2914      ! Allocate memory
2915      ALLOCATE( &
2916         & zts(     jpi, jpj, jpk, jpts ),  &
2917         & zts_adout( jpi, jpj, jpk, jpts ),  &
2918         & zrd_adin( jpi, jpj, jpk ),  &
2919         & zrd_adout(jpi, jpj, jpk ),  &
2920         & zts_tlin(  jpi, jpj, jpk, jpts ),  &
2921         & znts(      jpi, jpj, jpk, jpts ),  &
2922         & zts_tlout( jpi, jpj, jpk, jpts ),  &
2923         & zrd_tlout(jpi, jpj, jpk )    )
2924
2925      ! Initialize random field standard deviationsthe reference state
2926      zts = tsn
2927      ! store initial nn_eos
2928      jeos = nn_eos
2929      DO jn = 0, 2
2930         nn_eos = jn
2931         !=============================================================
2932         ! 1) dx = ( T ) and dy = ( T )
2933         !=============================================================
2934
2935         !--------------------------------------------------------------------
2936         ! Reset the tangent and adjoint variables
2937         !--------------------------------------------------------------------
2938         zts_tlin(:,:,:,:) = 0.0_wp
2939         zts_tlout(:,:,:,:) = 0.0_wp
2940         zrd_tlout(:,:,:) = 0.0_wp
2941         zts_adout(:,:,:,:) = 0.0_wp
2942         zrd_adin(:,:,:) = 0.0_wp
2943         zrd_adout(:,:,:) = 0.0_wp
2944
2945         !--------------------------------------------------------------------
2946         ! Initialize the tangent input with random noise: dx
2947         !--------------------------------------------------------------------
2948         CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2949         CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2950         DO jk = 1, jpk
2951            DO jj = nldj, nlej
2952               DO ji = nldi, nlei
2953                  zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2954               END DO
2955            END DO
2956         END DO
2957         !--------------------------------------------------------------------
2958         ! Call the tangent routine: dy = L dx
2959         !--------------------------------------------------------------------
2960         zts_tlout(:,:,:,:) = zts_tlin
2961         CALL eos_bn2_tan( zts, zts_tlout, zrd_tlout )
2962         !--------------------------------------------------------------------
2963         ! Initialize the adjoint variables: dy^* = W dy
2964         !--------------------------------------------------------------------
2965         DO jk = 1, jpk
2966            DO jj = nldj, nlej
2967               DO ji = nldi, nlei
2968                  zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2969                       &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
2970                       &               * tmask(ji,jj,jk)
2971               END DO
2972            END DO
2973         END DO
2974
2975         !--------------------------------------------------------------------
2976         ! Compute the scalar product: ( L dx )^T W dy
2977         !--------------------------------------------------------------------
2978
2979         zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2980
2981         !--------------------------------------------------------------------
2982         ! Call the adjoint routine: dx^* = L^T dy^*
2983         !--------------------------------------------------------------------
2984
2985         zrd_adout(:,:,:) = zrd_adin(:,:,:)
2986         CALL eos_bn2_adj( zts, zts_adout, zrd_adout )
2987
2988         !--------------------------------------------------------------------
2989         ! Compute the scalar product: dx^T L^T W dy
2990         !--------------------------------------------------------------------
2991         zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem )) + &
2992              & DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal ))
2993
2994         ! Compare the scalar products
2995         ! 14 char:'12345678901234'
2996         SELECT CASE( jn )
2997         CASE (0) ; cl_name = 'bn2_adj     T1'
2998         CASE (1) ; cl_name = 'bn2_adj     T2'
2999         CASE (2) ; cl_name = 'bn2_adj     T3'
3000         END SELECT
3001         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
3002      ENDDO
3003      ! restore initial nn_eos
3004      nn_eos = jeos
3005
3006      ! Deallocate memory
3007      DEALLOCATE(      &
3008         & zts,       &
3009         & zts_adout,   &
3010         & zrd_adin,   &
3011         & zrd_adout,  &
3012         & zts_tlin,    &
3013         & zts_tlout,   &
3014         & zrd_tlout,  &
3015         & znts    )
3016   END SUBROUTINE bn2_adj_tst
3017#endif
3018   !!======================================================================
3019END MODULE eosbn2_tam
Note: See TracBrowser for help on using the repository browser.