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 2840 for branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 – NEMO

Ignore:
Timestamp:
2011-09-15T14:38:15+02:00 (13 years ago)
Author:
agn
Message:

branches/2011/dev_r2782_NOCS_Griffies ticket #838 bugfixes + preliminary support for s-coords

File:
1 edited

Legend:

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

    r2715 r2840  
    33   !!                       ***  MODULE  eosbn2  *** 
    44   !! Ocean diagnostic variable : equation of state - in situ and potential density 
    5    !!                                               - Brunt-Vaisala frequency  
     5   !!                                               - Brunt-Vaisala frequency 
    66   !!============================================================================== 
    77   !! History :  OPA  ! 1989-03  (O. Marti)  Original code 
     
    2727   !!   eos_insitu_2d  : Compute the in situ density for 2d fields 
    2828   !!   eos_bn2        : Compute the Brunt-Vaisala frequency 
    29    !!   eos_alpbet     : calculates the in situ thermal and haline expansion coeff. 
     29   !!   eos_alpbet     : calculates the in situ thermal/haline expansion ratio 
    3030   !!   tfreez         : Compute the surface freezing temperature 
    3131   !!   eos_init       : set eos parameters (namelist) 
     
    4141   PRIVATE 
    4242 
    43    !                   !! * Interface  
     43   !                   !! * Interface 
    4444   INTERFACE eos 
    4545      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 
    46    END INTERFACE  
     46   END INTERFACE 
    4747   INTERFACE bn2 
    4848      MODULE PROCEDURE eos_bn2 
    49    END INTERFACE  
     49   END INTERFACE 
    5050 
    5151   PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules 
     
    6161 
    6262   REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio 
    63     
     63 
    6464   !! * Substitutions 
    6565#  include "domzgr_substitute.h90" 
     
    7575      !!---------------------------------------------------------------------- 
    7676      !!                   ***  ROUTINE eos_insitu  *** 
    77       !!  
    78       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from  
     77      !! 
     78      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    7979      !!       potential temperature and salinity using an equation of state 
    8080      !!       defined through the namelist parameter nn_eos. 
     
    134134!CDIR NOVERRCHK 
    135135         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    136          !   
     136         ! 
    137137         DO jk = 1, jpkm1 
    138138            DO jj = 1, jpj 
     
    199199      !!---------------------------------------------------------------------- 
    200200      !!                  ***  ROUTINE eos_insitu_pot  *** 
    201       !!            
     201      !! 
    202202      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the 
    203203      !!      potential volumic mass (Kg/m3) from potential temperature and 
    204       !!      salinity fields using an equation of state defined through the  
     204      !!      salinity fields using an equation of state defined through the 
    205205      !!     namelist parameter nn_eos. 
    206206      !! 
     
    230230      !!      nn_eos = 2 : linear equation of state function of temperature and 
    231231      !!               salinity 
    232       !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0  
     232      !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0 
    233233      !!                       = rn_beta * s - rn_alpha * tn - 1. 
    234234      !!              rhop(t,s)  = rho(t,s) 
     
    265265!CDIR NOVERRCHK 
    266266         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    267          !   
     267         ! 
    268268         DO jk = 1, jpkm1 
    269269            DO jj = 1, jpj 
     
    336336      !!                  ***  ROUTINE eos_insitu_2d  *** 
    337337      !! 
    338       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from  
     338      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    339339      !!      potential temperature and salinity using an equation of state 
    340340      !!      defined through the namelist parameter nn_eos. * 2D field case 
     
    374374      !                                                           ! 2 : salinity               [psu] 
    375375      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m] 
    376       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density  
     376      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
    377377      !! 
    378378      INTEGER  ::   ji, jj                    ! dummy loop indices 
     
    449449         DO jj = 1, jpjm1 
    450450            DO ji = 1, fs_jpim1   ! vector opt. 
    451                prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)  
     451               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
    452452            END DO 
    453453         END DO 
     
    468468      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the time- 
    469469      !!      step of the input arguments 
    470       !!       
     470      !! 
    471471      !! ** Method : 
    472472      !!       * nn_eos = 0  : UNESCO sea water properties 
     
    482482      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 
    483483      !!      The use of potential density to compute N^2 introduces e r r o r 
    484       !!      in the sign of N^2 at great depths. We recommand the use of  
     484      !!      in the sign of N^2 at great depths. We recommand the use of 
    485485      !!      nn_eos = 0, except for academical studies. 
    486486      !!        Macro-tasked on horizontal slab (jk-loop) 
     
    497497      !! 
    498498      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    499       REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars  
     499      REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars 
    500500#if defined key_zdfddm 
    501501      REAL(wp) ::   zds   ! local scalars 
     
    504504 
    505505      ! pn2 : interior points only (2=< jk =< jpkm1 ) 
    506       ! --------------------------  
     506      ! -------------------------- 
    507507      ! 
    508508      SELECT CASE( nn_eos ) 
     
    542542                     &                                - 0.121555e-07_wp ) * zh 
    543543                     ! 
    544                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2  
     544                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
    545545                     &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    546546                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
     
    565565               &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   & 
    566566               &               / fse3w(:,:,jk) * tmask(:,:,jk) 
    567          END DO  
     567         END DO 
    568568#if defined key_zdfddm 
    569569         DO jk = 2, jpkm1                                 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    570570            DO jj = 1, jpj 
    571571               DO ji = 1, jpi 
    572                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )   
     572                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 
    573573                  IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 
    574574                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     
    587587 
    588588 
    589    SUBROUTINE eos_alpbet( pts, palph, pbeta ) 
    590       !!---------------------------------------------------------------------- 
    591       !!                 ***  ROUTINE ldf_slp_grif  *** 
    592       !! 
    593       !! ** Purpose :   Calculates the thermal and haline expansion coefficients at T-points.  
    594       !! 
    595       !! ** Method  :   calculates alpha and beta at T-points  
     589   SUBROUTINE eos_alpbet( pts, palpbet, beta0 ) 
     590      !!---------------------------------------------------------------------- 
     591      !!                 ***  ROUTINE eos_alpbet  *** 
     592      !! 
     593      !! ** Purpose :   Calculates the in situ thermal/haline expansion ratio at T-points 
     594      !! 
     595      !! ** Method  :   calculates alpha / beta ratio at T-points 
    596596      !!       * nn_eos = 0  : UNESCO sea water properties 
    597       !!         The brunt-vaisala frequency is computed using the polynomial 
    598       !!      polynomial expression of McDougall (1987): 
    599       !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 
    600       !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
    601       !!      computed and used in zdfddm module : 
    602       !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 
     597      !!                       The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 
     598      !!                       polynomial expression of McDougall (1987). 
     599      !!                       Scalar beta0 is returned = 1. 
    603600      !!       * nn_eos = 1  : linear equation of state (temperature only) 
    604       !!            N^2 = grav * rn_alpha * dk[ t ]/e3w 
     601      !!                       The ratio is undefined, so we return alpha as palpbet 
     602      !!                       Scalar beta0 is returned = 0. 
    605603      !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    606       !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 
    607       !!       * nn_eos = 3  : Jackett JAOT 2003 ??? 
    608       !! 
    609       !! ** Action  : - palph, pbeta : thermal and haline expansion coeff. at T-point 
    610       !!---------------------------------------------------------------------- 
    611       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts            ! pot. temperature & salinity 
    612       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palph, pbeta   ! thermal & haline expansion coeff. 
    613       ! 
     604      !!                       The alpha/beta ratio is returned as ralpbet 
     605      !!                       Scalar beta0 is returned = 1. 
     606      !! 
     607      !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points 
     608      !!            :   beta0   : 1. or 0. 
     609      !!---------------------------------------------------------------------- 
     610      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity 
     611      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpbet   ! thermal/haline expansion ratio 
     612      REAL(wp),                              INTENT(  out) ::   beta0     ! set = 1 except with case 1 eos, rho=rho(T) 
     613      !! 
    614614      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    615       REAL(wp) ::   zt, zs, zh   ! local scalars  
     615      REAL(wp) ::   zt, zs, zh   ! local scalars 
    616616      !!---------------------------------------------------------------------- 
    617617      ! 
     
    624624                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature 
    625625                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
    626                   zh = fsdept(ji,jj,jk)              ! depth in meters  
    627                   ! 
    628                   pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt   & 
    629                      &                                        - 0.301985e-05_wp ) * zt   & 
    630                      &           + 0.785567e-03_wp                                       & 
    631                      &           + (     0.515032e-08_wp * zs                            & 
    632                      &                 + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs   & 
    633                      &           + ( (   0.121551e-17_wp * zh                            & 
    634                      &                 - 0.602281e-15_wp * zs                            & 
    635                      &                 - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh   & 
    636                      &                                        + 0.408195e-10_wp   * zs   & 
    637                      &             + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt   & 
    638                      &                                        - 0.121555e-07_wp ) * zh 
    639                      ! 
    640                   palph(ji,jj,jk) = - pbeta(ji,jj,jk) *                             & 
    641                       &     ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   & 
    642                       &                                  - 0.203814e-03_wp ) * zt   & 
    643                       &                                  + 0.170907e-01_wp ) * zt   & 
    644                       &   + 0.665157e-01_wp                                         & 
    645                       &   +     ( - 0.678662e-05_wp * zs                            & 
    646                       &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
    647                       &   +   ( ( - 0.302285e-13_wp * zh                            & 
    648                       &           - 0.251520e-11_wp * zs                            & 
    649                       &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
    650                       &           - 0.164759e-06_wp * zs                            & 
    651                       &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
    652                       &                                  + 0.380374e-04_wp ) * zh) 
     626                  zh = fsdept(ji,jj,jk)               ! depth in meters 
     627                  ! 
     628                  palpbet(ji,jj,jk) =                                              & 
     629                     &     ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   & 
     630                     &                                  - 0.203814e-03_wp ) * zt   & 
     631                     &                                  + 0.170907e-01_wp ) * zt   & 
     632                     &   + 0.665157e-01_wp                                         & 
     633                     &   +     ( - 0.678662e-05_wp * zs                            & 
     634                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
     635                     &   +   ( ( - 0.302285e-13_wp * zh                            & 
     636                     &           - 0.251520e-11_wp * zs                            & 
     637                     &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
     638                     &           - 0.164759e-06_wp * zs                            & 
     639                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
     640                     &                                  + 0.380374e-04_wp ) * zh 
    653641               END DO 
    654642            END DO 
    655643         END DO 
    656          ! 
    657       CASE ( 1 ) 
    658          palph(:,:,:) = - rn_alpha 
    659          pbeta(:,:,:) =   0._wp 
    660          ! 
    661       CASE ( 2 ) 
    662          palph(:,:,:) = - rn_alpha 
    663          pbeta(:,:,:) =   rn_beta 
     644         beta0 = 1._wp 
     645         ! 
     646      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==! 
     647         palpbet(:,:,:) = rn_alpha 
     648         beta0 = 0._wp 
     649         ! 
     650      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
     651         palpbet(:,:,:) = ralpbet 
     652         beta0 = 1._wp 
    664653         ! 
    665654      CASE DEFAULT 
     
    747736 
    748737   !!====================================================================== 
    749 END MODULE eosbn2   
     738END MODULE eosbn2 
Note: See TracChangeset for help on using the changeset viewer.