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.
Changeset 2371 for branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 – NEMO

Ignore:
Timestamp:
2010-11-10T16:38:27+01:00 (13 years ago)
Author:
acc
Message:

nemo_v3_3_beta. Improvement of the Griffies operator code (#680). Some aspects are still undergoing scientific assessment but the code structure is ready for release. Dynamical allocation of the large triad arrays has been introduced to reduce memory when the new operator is not in use.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r2287 r2371  
    1717   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function 
    1818   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
     19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    2627   !!   eos_insitu_2d  : Compute the in situ density for 2d fields 
    2728   !!   eos_bn2        : Compute the Brunt-Vaisala frequency 
     29   !!   eos_alpbet     : calculates the in situ thermal and haline expansion coeff. 
    2830   !!   tfreez         : Compute the surface freezing temperature 
    2931   !!   eos_init       : set eos parameters (namelist) 
     
    3840   PRIVATE 
    3941 
    40    !! * Interface  
     42   !                   !! * Interface  
    4143   INTERFACE eos 
    4244      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 
     
    4951   PUBLIC   eos_init   ! called by istate module 
    5052   PUBLIC   bn2        ! called by step module 
     53   PUBLIC   eos_alpbet ! called by ldfslp module 
    5154   PUBLIC   tfreez     ! called by sbcice_... modules 
    5255 
    53    !                                        !!* Namelist (nameos) * 
    54    INTEGER , PUBLIC ::   nn_eos   = 0        !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
    55    REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4  !: thermal expension coeff. (linear equation of state) 
    56    REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4  !: saline  expension coeff. (linear equation of state) 
    57  
    58    REAL(wp), PUBLIC ::   ralpbet           !: alpha / beta ratio 
     56   !                                          !!* Namelist (nameos) * 
     57   INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
     58   REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state) 
     59   REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4_wp !: saline  expension coeff. (linear equation of state) 
     60 
     61   REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio 
    5962    
    6063   !! * Substitutions 
     
    6467   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    6568   !! $Id$ 
    66    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     69   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6770   !!---------------------------------------------------------------------- 
    68  
    6971CONTAINS 
    7072 
     
    135137                  ! 
    136138                  ! compute volumic mass pure water at atm pressure 
    137                   zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt   & 
    138                      &      -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 
     139                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     140                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    139141                  ! seawater volumic mass atm pressure 
    140                   zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt        & 
    141                      &                   -4.0899e-3 ) *zt+0.824493 
    142                   zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 
    143                   zr4= 4.8314e-4 
     142                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     143                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     144                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     145                  zr4= 4.8314e-4_wp 
    144146                  ! 
    145147                  ! potential volumic mass (reference to the surface) 
     
    147149                  ! 
    148150                  ! add the compression terms 
    149                   ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 
    150                   zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 
     151                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     152                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    151153                  zb = zbw + ze * zs 
    152154                  ! 
    153                   zd = -2.042967e-2 
    154                   zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 
    155                   zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 
     155                  zd = -2.042967e-2_wp 
     156                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
     157                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    156158                  za = ( zd*zsr + zc ) *zs + zaw 
    157159                  ! 
    158                   zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545 
    159                   za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 
    160                   zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 
     160                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     161                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     162                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    161163                  zk0= ( zb1*zsr + za1 )*zs + zkw 
    162164                  ! 
    163165                  ! masked in situ density anomaly 
    164                   prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
     166                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    165167                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk) 
    166168               END DO 
     
    170172      CASE( 1 )                !==  Linear formulation function of temperature only  ==! 
    171173         DO jk = 1, jpkm1 
    172             prd(:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
     174            prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
    173175         END DO 
    174176         ! 
     
    258260                  ! 
    259261                  ! compute volumic mass pure water at atm pressure 
    260                   zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt   & 
    261                      &                       -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 
     262                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
     263                     &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
    262264                  ! seawater volumic mass atm pressure 
    263                   zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt   & 
    264                      &                                   -4.0899e-3 ) *zt+0.824493 
    265                   zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 
    266                   zr4= 4.8314e-4 
     265                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   & 
     266                     &                                         -4.0899e-3_wp ) *zt+0.824493_wp 
     267                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     268                  zr4= 4.8314e-4_wp 
    267269                  ! 
    268270                  ! potential volumic mass (reference to the surface) 
     
    273275                  ! 
    274276                  ! add the compression terms 
    275                   ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 
    276                   zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 
     277                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     278                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    277279                  zb = zbw + ze * zs 
    278280                  ! 
    279                   zd = -2.042967e-2 
    280                   zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 
    281                   zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 
     281                  zd = -2.042967e-2_wp 
     282                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
     283                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
    282284                  za = ( zd*zsr + zc ) *zs + zaw 
    283285                  ! 
    284                   zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545 
    285                   za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 
    286                   zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 
     286                  zb1=   (  -0.1909078_wp  *zt+7.390729_wp    ) *zt-55.87545_wp 
     287                  za1= ( (   2.326469e-3_wp*zt+1.553190_wp    ) *zt-65.00517_wp ) *zt + 1044.077_wp 
     288                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    287289                  zk0= ( zb1*zsr + za1 )*zs + zkw 
    288290                  ! 
    289291                  ! masked in situ density anomaly 
    290                   prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
     292                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    291293                     &             - rau0  ) * zrau0r * tmask(ji,jj,jk) 
    292294               END DO 
     
    296298      CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    297299         DO jk = 1, jpkm1 
    298             prd  (:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_sal) )        * tmask(:,:,jk) 
    299             prhop(:,:,jk) = ( 1.e0   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk) 
     300            prd  (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_sal) )        * tmask(:,:,jk) 
     301            prhop(:,:,jk) = ( 1.e0_wp   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk) 
    300302         END DO 
    301303         ! 
     
    303305         DO jk = 1, jpkm1 
    304306            prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk) 
    305             prhop(:,:,jk) = ( 1.e0     + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk) 
     307            prhop(:,:,jk) = ( 1.e0_wp  + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk) 
    306308         END DO 
    307309         ! 
     
    382384               ! 
    383385               ! compute volumic mass pure water at atm pressure 
    384                zr1 = ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt   & 
    385                   &                        -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 
     386               zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
     387                  &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
    386388               ! seawater volumic mass atm pressure 
    387                zr2 = ( ( ( 5.3875e-9*zt-8.2467e-7 )*zt+7.6438e-5 ) *zt   & 
    388                   &                                   -4.0899e-3 ) *zt+0.824493 
    389                zr3 = ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 
    390                zr4 = 4.8314e-4 
     389               zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   & 
     390                  &                                   -4.0899e-3_wp ) *zt+0.824493_wp 
     391               zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 
     392               zr4 = 4.8314e-4_wp 
    391393               ! 
    392394               ! potential volumic mass (reference to the surface) 
     
    394396               ! 
    395397               ! add the compression terms 
    396                ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 
    397                zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 
     398               ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     399               zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    398400               zb = zbw + ze * zs 
    399401               ! 
    400                zd = -2.042967e-2 
    401                zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 
    402                zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt -4.721788 
     402               zd =    -2.042967e-2_wp 
     403               zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
     404               zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp 
    403405               za = ( zd*zsr + zc ) *zs + zaw 
    404406               ! 
    405                zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545 
    406                za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 
    407                zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt   & 
    408                   &                                       +2098.925 ) *zt+190925.6 
     407               zb1=     (-0.1909078_wp  *zt+7.390729_wp      ) *zt-55.87545_wp 
     408               za1=   ( ( 2.326469e-3_wp*zt+1.553190_wp      ) *zt-65.00517_wp ) *zt+1044.077_wp 
     409               zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp   ) *zt-30.41638_wp ) *zt   & 
     410                  &                             +2098.925_wp ) *zt+190925.6_wp 
    409411               zk0= ( zb1*zsr + za1 )*zs + zkw 
    410412               ! 
    411413               ! masked in situ density anomaly 
    412                prd(ji,jj) = ( zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
     414               prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
    413415            END DO 
    414416         END DO 
     
    417419         DO jj = 1, jpjm1 
    418420            DO ji = 1, fs_jpim1   ! vector opt. 
    419                prd(ji,jj) = ( 0.0285 - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
     421               prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
    420422            END DO 
    421423         END DO 
     
    490492                  zh = fsdepw(ji,jj,jk)                                     ! depth in meters  at w-point 
    491493                  ! 
    492                   zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta 
    493                      &                               - 0.203814e-03 ) * zt   & 
    494                      &                               + 0.170907e-01 ) * zt   & 
    495                      &   + 0.665157e-01                                      & 
    496                      &   +     ( - 0.678662e-05 * zs                         & 
    497                      &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   & 
    498                      &   +   ( ( - 0.302285e-13 * zh                         & 
    499                      &           - 0.251520e-11 * zs                         & 
    500                      &           + 0.512857e-12 * zt * zt           ) * zh   & 
    501                      &           - 0.164759e-06 * zs                         & 
    502                      &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   & 
    503                      &                               + 0.380374e-04 ) * zh 
     494                  zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &   ! ratio alpha/beta 
     495                     &                                  - 0.203814e-03_wp ) * zt   & 
     496                     &                                  + 0.170907e-01_wp ) * zt   & 
     497                     &   +         0.665157e-01_wp                                 & 
     498                     &   +     ( - 0.678662e-05_wp * zs                            & 
     499                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
     500                     &   +   ( ( - 0.302285e-13_wp * zh                            & 
     501                     &           - 0.251520e-11_wp * zs                            & 
     502                     &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
     503                     &           - 0.164759e-06_wp * zs                            & 
     504                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
     505                     &                                  + 0.380374e-04_wp ) * zh 
    504506                     ! 
    505                   zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta 
    506                      &                            - 0.301985e-05 ) * zt      & 
    507                      &   + 0.785567e-03                                      & 
    508                      &   + (     0.515032e-08 * zs                           & 
    509                      &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     & 
    510                      &   +(  (   0.121551e-17 * zh                           & 
    511                      &         - 0.602281e-15 * zs                           & 
    512                      &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     & 
    513                      &                             + 0.408195e-10   * zs     & 
    514                      &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     & 
    515                      &                             - 0.121555e-07 ) * zh 
     507                  zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta 
     508                     &                               - 0.301985e-05_wp ) * zt      & 
     509                     &   +       0.785567e-03_wp                                   & 
     510                     &   + (     0.515032e-08_wp * zs                              & 
     511                     &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     & 
     512                     &   + ( (   0.121551e-17_wp * zh                              & 
     513                     &         - 0.602281e-15_wp * zs                              & 
     514                     &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     & 
     515                     &                                + 0.408195e-10_wp   * zs     & 
     516                     &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     & 
     517                     &                                - 0.121555e-07_wp ) * zh 
    516518                     ! 
    517519                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2  
     
    521523                  !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!! 
    522524                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    523                   IF ( ABS( zds) <= 1.e-20 ) zds = 1.e-20 
     525                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    524526                  rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    525527#endif 
     
    544546               DO ji = 1, jpi 
    545547                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )   
    546                   IF ( ABS( zds ) <= 1.e-20 ) zds = 1.e-20 
     548                  IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 
    547549                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    548550               END DO 
     
    560562 
    561563 
     564   SUBROUTINE eos_alpbet( pts, palph, pbeta ) 
     565      !!---------------------------------------------------------------------- 
     566      !!                 ***  ROUTINE ldf_slp_grif  *** 
     567      !! 
     568      !! ** Purpose :   Calculates the thermal and haline expansion coefficients at T-points.  
     569      !! 
     570      !! ** Method  :   calculates alpha and beta at T-points  
     571      !!       * nn_eos = 0  : UNESCO sea water properties 
     572      !!         The brunt-vaisala frequency is computed using the polynomial 
     573      !!      polynomial expression of McDougall (1987): 
     574      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 
     575      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
     576      !!      computed and used in zdfddm module : 
     577      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 
     578      !!       * nn_eos = 1  : linear equation of state (temperature only) 
     579      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w 
     580      !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
     581      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 
     582      !!       * nn_eos = 3  : Jackett JAOT 2003 ??? 
     583      !! 
     584      !! ** Action  : - palph, pbeta : thermal and haline expansion coeff. at T-point 
     585      !!---------------------------------------------------------------------- 
     586      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts            ! pot. temperature & salinity 
     587      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palph, pbeta   ! thermal & haline expansion coeff. 
     588      !! 
     589      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     590      REAL(wp) ::   zt, zs, zh   ! local scalars  
     591      !!---------------------------------------------------------------------- 
     592      ! 
     593      SELECT CASE ( nn_eos ) 
     594      ! 
     595      CASE ( 0 )               ! Jackett and McDougall (1994) formulation 
     596         DO jk = 1, jpk 
     597            DO jj = 1, jpj 
     598               DO ji = 1, jpi 
     599                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature 
     600                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
     601                  zh = fsdept(ji,jj,jk)              ! depth in meters  
     602                  ! 
     603                  pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt   & 
     604                     &                                        - 0.301985e-05_wp ) * zt   & 
     605                     &           + 0.785567e-03_wp                                       & 
     606                     &           + (     0.515032e-08_wp * zs                            & 
     607                     &                 + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs   & 
     608                     &           + ( (   0.121551e-17_wp * zh                            & 
     609                     &                 - 0.602281e-15_wp * zs                            & 
     610                     &                 - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh   & 
     611                     &                                        + 0.408195e-10_wp   * zs   & 
     612                     &             + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt   & 
     613                     &                                        - 0.121555e-07_wp ) * zh 
     614                     ! 
     615                  palph(ji,jj,jk) = - pbeta(ji,jj,jk) *                             & 
     616                      &     ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   & 
     617                      &                                  - 0.203814e-03_wp ) * zt   & 
     618                      &                                  + 0.170907e-01_wp ) * zt   & 
     619                      &   + 0.665157e-01_wp                                         & 
     620                      &   +     ( - 0.678662e-05_wp * zs                            & 
     621                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
     622                      &   +   ( ( - 0.302285e-13_wp * zh                            & 
     623                      &           - 0.251520e-11_wp * zs                            & 
     624                      &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
     625                      &           - 0.164759e-06_wp * zs                            & 
     626                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
     627                      &                                  + 0.380374e-04_wp ) * zh) 
     628               END DO 
     629            END DO 
     630         END DO 
     631         ! 
     632      CASE ( 1 ) 
     633         palph(:,:,:) = - rn_alpha 
     634         pbeta(:,:,:) =   0._wp 
     635         ! 
     636      CASE ( 2 ) 
     637         palph(:,:,:) = - rn_alpha 
     638         pbeta(:,:,:) =   rn_beta 
     639         ! 
     640      CASE DEFAULT 
     641         IF(lwp) WRITE(numout,cform_err) 
     642         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     643         nstop = nstop + 1 
     644         ! 
     645      END SELECT 
     646      ! 
     647   END SUBROUTINE eos_alpbet 
     648 
     649 
    562650   FUNCTION tfreez( psal ) RESULT( ptf ) 
    563651      !!---------------------------------------------------------------------- 
     
    576664      !!---------------------------------------------------------------------- 
    577665      ! 
    578       ptf(:,:) = ( - 0.0575 + 1.710523e-3 * SQRT( psal(:,:) )   & 
    579          &                  - 2.154996e-4 *       psal(:,:)   ) * psal(:,:) 
     666      ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
     667         &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:) 
    580668      ! 
    581669   END FUNCTION tfreez 
Note: See TracChangeset for help on using the changeset viewer.