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 3894 – NEMO

Changeset 3894


Ignore:
Timestamp:
2013-04-27T17:46:22+02:00 (11 years ago)
Author:
gm
Message:

dev_r3858_CNRS3_Ediag: #927 simplification of eosbn2, use of alpha & beta in zdfddm, trabbl...

Location:
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM
Files:
24 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

    r3893 r3894  
    877877   ln_glo_trd  = .FALSE.   ! (T) global domain averaged diag for T, T^2, KE, and PE 
    878878   ln_dyn_trd  = .FALSE.   ! (T) 3D momentum trend output 
    879    ln_dyn_mld  = .FALSE.   ! (T) 2D momentum trends averaged over the mixed layer  
     879   ln_dyn_mxl  = .FALSE.   ! (T) 2D momentum trends averaged over the mixed layer  
    880880   ln_vor_trd  = .FALSE.   ! (T) 2D barotropic vorticity trends  
    881881   ln_KE_trd   = .FALSE.   ! (T) 3D Kinetic   Energy     trends 
    882882   ln_PE_trd   = .FALSE.   ! (T) 3D Potential Energy     trends 
    883883   ln_tra_trd  = .FALSE.   ! (T) 3D tracer trend output 
    884    ln_tra_mld  = .FALSE.   ! (T) 2D tracer trends averaged over the mixed layer  
     884   ln_tra_mxl  = .FALSE.   ! (T) 2D tracer trends averaged over the mixed layer  
    885885   nn_trd      = 365       !  print frequency (ln_glo_trd=T) (unit=time step) 
    886886/ 
     
    894894   cn_trdrst_in      = "restart_mld"   ! suffix of ocean restart name (input) 
    895895   cn_trdrst_out     = "restart_mld"   ! suffix of ocean restart name (output) 
    896    ln_trdmld_restart = .false.         !  restart for ML diagnostics 
    897    ln_trdmld_instant = .false.         !  flag to diagnose trends of instantantaneous or mean ML T/S 
     896   ln_trdmxl_restart = .false.         !  restart for ML diagnostics 
     897   ln_trdmxl_instant = .false.         !  flag to diagnose trends of instantantaneous or mean ML T/S 
    898898/ 
    899899!----------------------------------------------------------------------- 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90

    r3625 r3894  
    7171      IF( .NOT. ln_limini ) THEN   
    7272          
    73          tfu(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
     73         tfu(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
    7474 
    7575         DO jj = 1, jpj 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r3625 r3894  
    9494      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 
    9595 
    96       t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
     96      t_bo(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
    9797 
    9898      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r3827 r3894  
    243243      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    244244      ! 
    245       CALL eos    ( tsn, rhd, rhop )                                       ! In any case, we need rhop 
    246       CALL zdf_mxl( kt )                                                   ! In any case, we need mxl  
     245      CALL eos( tsn, fsdept(:,:,:), rhd, rhop )                         ! In any case, we need rhop 
     246      CALL zdf_mxl( kt )                                                ! In any case, we need mxl  
    247247      ! 
    248248      avt(:,:,:)       = sf_dyn(jf_avt)%fnow(:,:,:)  * tmask(:,:,:)    ! vertical diffusive coefficient  
     
    546546      !!--------------------------------------------------------------------- 
    547547#if defined key_ldfslp && ! defined key_c1d 
    548       CALL eos( pts, rhd, rhop )   ! Time-filtered in situ density  
    549       CALL bn2( pts, rn2 )         ! before Brunt-Vaisala frequency 
     548      CALL eos( pts, fsdept(:,:,:), rhd, rhop )   ! Time-filtered in situ density  
     549      CALL eos( tsn, rab_n, rn2  )                ! before Brunt-Vaisala frequency 
    550550      IF( ln_zps )   & 
    551551         &  CALL zps_hde( kt, jpts, pts, gtsu, gtsv, rhd, gru, grv )  ! Partial steps: before Horizontal DErivative 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r3785 r3894  
    679679      !!---------------------------------------------------------------------- 
    680680 
    681       ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)  
     681      ! freezing point 
    682682      ! used to prevent the applied increments taking the temperature below the local freezing point  
    683  
    684683#if defined key_cice  
    685684        fzptnz(:,:,:) = -1.8_wp 
    686685#else  
    687         DO jk = 1, jpk 
    688            DO jj = 1, jpj 
    689               DO ji = 1, jpk 
    690                  fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) )                   &  
    691                                                   - 2.154996e-4_wp *       tsn(ji,jj,jk,jp_sal)   ) * tsn(ji,jj,jk,jp_sal)  &  
    692                                                   - 7.53e-4_wp * fsdepw(ji,jj,jk)       ! (pressure in dbar)  
    693               END DO 
    694            END DO 
    695         END DO 
     686      DO jk = 1, jpk 
     687         fzptnz(:,:,jk) =eos_fzp( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) )  
     688      END DO 
    696689#endif  
    697690 
     
    709702            IF(lwp) THEN 
    710703               WRITE(numout,*)  
    711                WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', & 
    712                   &  kt,' with IAU weight = ', wgtiau(it) 
     704               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 
    713705               WRITE(numout,*) '~~~~~~~~~~~~' 
    714706            ENDIF 
     
    758750            IF (ln_temnofreeze) THEN 
    759751               ! Do not apply negative increments if the temperature will fall below freezing 
    760                WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 
    761                   &   tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )  
     752               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )  
    762753                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
    763754               END WHERE 
     
    768759               ! Do not apply negative increments if the salinity will fall below a specified 
    769760               ! minimum value salfixmin 
    770                WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 
    771                   &   tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )  
     761               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )  
    772762                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
    773763               END WHERE 
     
    778768            tsb(:,:,:,:) = tsn(:,:,:,:)               ! Update before fields 
    779769 
    780             CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities 
     770            CALL eos( tsb, fsdept(:,:,:), rhd, rhop )                ! Before potential and in situ densities 
    781771          
    782772            IF( ln_zps .AND. .NOT. lk_c1d ) & 
     
    786776 
    787777#if defined key_zdfkpp 
    788             CALL eos( tsn, rhd )                      ! Compute rhd 
     778            CALL eos( tsn, fsdept(:,:,:), rhd )                      ! Compute rhd 
    789779#endif 
    790780 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90

    r3680 r3894  
    7070      ! Ocean physics update                (ua, va, ta, sa used as workspace) 
    7171      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    72                          CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency 
    73                          CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency 
     72      !                                               ! thermal/haline expansion coef. & Brunt-Vaisala frequency 
     73                         CALL eos( tsb, rab_b, rn2b )    ! before  
     74                         CALL eos( tsn, rab_n, rn2  )    ! now 
     75      ! 
    7476      !  VERTICAL PHYSICS    
    7577                         CALL zdf_bfr( kstp )         ! bottom friction 
     
    125127                             CALL tra_nxt    ( kstp )        ! tracer fields at next time step 
    126128      IF( ln_zdfnpc      )   CALL tra_npc    ( kstp )        ! applied non penetrative convective adjustment on (t,s) 
    127                              CALL eos( tsb, rhd, rhop )      ! now (swap=before) in situ density for dynhpg module 
     129                                                             ! now (swap=before) in situ density for dynhpg module 
     130                             CALL eos( tsb, fsdept(:,:,:), rhd, rhop ) 
    128131 
    129132      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    130133      ! Dynamics                                    (ta, sa used as workspace) 
    131134      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    132                                ua(:,:,:) = 0.e0               ! set dynamics trends to zero 
    133                                va(:,:,:) = 0.e0 
     135                             ua(:,:,:) = 0.e0               ! set dynamics trends to zero 
     136                             va(:,:,:) = 0.e0 
    134137 
    135                                CALL dyn_cor_c1d( kstp )       ! vorticity term including Coriolis 
    136                                CALL dyn_zdf    ( kstp )       ! vertical diffusion 
    137                                CALL dyn_nxt_c1d( kstp )       ! lateral velocity at next time step 
     138                             CALL dyn_cor_c1d( kstp )       ! vorticity term including Coriolis 
     139                             CALL dyn_zdf    ( kstp )       ! vertical diffusion 
     140                             CALL dyn_nxt_c1d( kstp )       ! lateral velocity at next time step 
    138141 
    139142      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    141144      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    142145                                 CALL stp_ctl( kstp, indic ) 
    143       IF( kstp == nit000     )   CALL iom_close( numror )             ! close input  ocean restart file 
    144       IF( lrst_oce           )   CALL rst_write  ( kstp )             ! write output ocean restart file 
     146      IF( kstp == nit000 )   CALL iom_close( numror )             ! close input  ocean restart file 
     147      IF( lrst_oce       )   CALL rst_write  ( kstp )             ! write output ocean restart file 
    145148      ! 
    146149   END SUBROUTINE stp_c1d 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r3294 r3894  
    9696 
    9797      !                      
    98       ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh 
     98      ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)    ! thermosteric ssh 
    9999      ztsn(:,:,:,jp_sal) = sn0(:,:,:) 
    100       CALL eos( ztsn, zrhd )                       ! now in situ density using initial salinity 
    101       ! 
    102       zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
     100      CALL eos( ztsn, fsdept(:,:,:), zrhd )         ! now in situ density using initial salinity 
     101      ! 
     102      zbotpres(:,:) = 0._wp                         ! no atmospheric surface pressure, levitating sea-ice 
    103103      DO jk = 1, jpkm1 
    104104         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 
     
    112112       
    113113      !                                         ! steric sea surface height 
    114       CALL eos( tsn, zrhd, zrhop )                 ! now in situ and potential density 
     114      CALL eos( tsn, fsdept(:,:,:), zrhd, zrhop )   ! now in situ and potential density 
    115115      zrhop(:,:,jpk) = 0._wp 
    116116      CALL iom_put( 'rhop', zrhop ) 
    117117      ! 
    118       zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
     118      zbotpres(:,:) = 0._wp                         ! no atmospheric surface pressure, levitating sea-ice 
    119119      DO jk = 1, jpkm1 
    120120         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r3893 r3894  
    122122         ENDIF 
    123123         ! 
    124          CALL eos( tsb, rhd, rhop )        ! before potential and in situ densities 
     124         CALL eos( tsb, fsdept(:,:,:), rhd, rhop )        ! before potential and in situ densities 
    125125#if ! defined key_c1d 
    126126         IF( ln_zps )   CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  & ! zps: before hor. gradient 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r3878 r3894  
    209209         CALL iom_get( numror, jpdom_autoglo, 'rhop'   , rhop    )   ! now    potential density 
    210210      ELSE 
    211          CALL eos    ( tsn, rhd, rhop )    
     211         CALL eos( tsn, fsdept(:,:,:), rhd, rhop )    
    212212      ENDIF 
    213213#if defined key_zdfkpp 
     
    215215         CALL iom_get( numror, jpdom_autoglo, 'rhd'    , rhd     )   ! now    in situ density anomaly 
    216216      ELSE 
    217          CALL eos( tsn, rhd )   ! compute rhd 
     217         CALL eos( tsn, fsdept(:,:,:), rhd )   ! compute rhd 
    218218      ENDIF 
    219219#endif 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90

    r3625 r3894  
    9898                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 ) 
    9999          
    100          fr_i(:,:) = tfreez( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius] 
     100         fr_i(:,:) = eos_fzp( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius] 
    101101#if defined key_coupled  
    102102         a_i(:,:,1) = fr_i(:,:)          
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r3625 r3894  
    122122         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean) 
    123123         ! 
    124          t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin] 
     124         t_bo(:,:) = eos_fzp( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin] 
    125125         !                                           ! (set to rt0 over land) 
    126126         CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90

    r3680 r3894  
    140140 
    141141         ! ... masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
    142          tfu(:,:) = tfreez( sss_m ) +  rt0  
     142         tfu(:,:) = eos_fzp( sss_m ) +  rt0  
    143143 
    144144         zsist (:,:,1) = sist (:,:) + rt0 * ( 1. - tmask(:,:,1) ) 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r3893 r3894  
    1515   !!             -   ! 2002-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d 
    1616   !!             -   ! 2003-08  (G. Madec)  F90, free form 
    17    !!            3.0  ! 2006-08  (G. Madec)  add tfreez function 
     17   !!            3.0  ! 2006-08  (G. Madec)  add tfreez function (now eos_fzp function) 
    1818   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
    19    !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_ab used in ldfslp 
    20    !!            3.5  ! 2012-03  (F. Roquet, G. Madec)  add pen_ddt_dds used in trdpen 
     19   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add alpha/beta used in ldfslp 
     20   !!            3.5  ! 2012-03  (F. Roquet, G. Madec)  add primitive of alpha and beta used in PE computation 
    2121   !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state 
    22    !!             -   ! 2012-05  (F. Roquet)  add eos_alpha_beta to be used in dynhpg 
     22   !!             -   ! 2013-04  (F. Roquet, G. Madec)  add eos_rab, change bn2 computation and reorganize the module 
    2323   !!---------------------------------------------------------------------- 
    2424 
    2525   !!---------------------------------------------------------------------- 
    2626   !!   eos            : generic interface of the equation of state 
    27    !!   eos_insitu     : density anomaly (rho/rho0 - 1) 
    28    !!   eos_insitu_pot : density and surface referenced potential density anomalies 
    29    !!   eos_insitu_2d  : density anomaly for 2d fields 
    30    !!   eos_bn2        : Brunt-Vaisala frequency 
    31    !!   eos_rab        : partial derivative of density anomaly with respect to T and S 
    32    !!   tfreez         : surface freezing temperature 
     27   !!      eos_rhd_pot : density anomaly (rho/rho0 - 1) and surface referenced potential density (T-point) 
     28   !!      eos_rab_bn2 : partial derivative of density anomaly with respect to T and S (T-point) 
     29   !!                    and Brunt-Vaisala frequency (w-point) 
     30   !!   eos_bn2        : Brunt-Vaisala frequency (w-point) computed from alpha & beta 
     31   !!   eos_fzp        : freezing temperature 
     32   !!   eos_pen        : primitive of partial derivative of density anomaly with respect to T and S (T-point) 
    3333   !!   eos_init       : set eos parameters (namelist) 
    3434   !!---------------------------------------------------------------------- 
     35   USE oce             ! ocean variables 
    3536   USE dom_oce         ! ocean space and time domain 
    3637   USE phycst          ! physical constants 
     38   USE ldfslp          ! lateral  physics: iso-neutral slopes 
    3739   USE zdfddm          ! vertical physics: double diffusion 
     40   ! 
    3841   USE in_out_manager  ! I/O manager 
    3942   USE lib_mpp         ! MPP library 
     
    4144   USE wrk_nemo        ! Memory Allocation 
    4245   USE timing          ! Timing 
     46   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    4347 
    4448   IMPLICIT NONE 
    4549   PRIVATE 
    4650 
    47    !                   !! * Interface 
    4851   INTERFACE eos 
    49       MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, & 
    50          &                eos_ab, eos_insitu_ab 
     52      MODULE PROCEDURE eos_rhd_pot, eos_rab_bn2 
    5153   END INTERFACE 
    52    INTERFACE bn2 
    53       MODULE PROCEDURE eos_bn2, eos_bn2_ab 
    54    END INTERFACE 
    55  
    56    PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules 
    57    PUBLIC   eos_init       ! called by istate module 
    58    PUBLIC   bn2            ! called by step module, update alpha/beta and compute bn2 
    59    PUBLIC   eos_pen        ! used for pe diagnostics in trdpen 
    60    PUBLIC   tfreez         ! called by sbcice_... modules 
    61  
    62    !                                          !!* Namelist (nameos) * 
    63    INTEGER , PUBLIC ::   nn_eos   = 0         !: = -1/0/1 type of eq. of state and Brunt-Vaisala frequ. 
    64  
    65    REAL(wp) ::   rn_alph0  = 1.67e-4 !: thermal expansion coeff. (simplified equation of state) 
    66    REAL(wp) ::   rn_beta0  = 7.8e-4 !: saline  expansion coeff. (simplified equation of state) 
    67    REAL(wp) ::   rn_gamm0  = 4.4e-6 !: depth expansion coeff. (simplified equation of state) 
    68    REAL(wp) ::   rn_cabbel = 3.0e-2 !: cabbeling coeff. (simplified equation of state) 
    69    REAL(wp) ::   rn_thermo = 1.1e-4 !: thermobaric coeff. (simplified equation of state) 
    70    REAL(wp) ::   offset    = 0      !: offset applied on Vallis density anomaly, so that rho(10,35,0)=1027 
     54 
     55   PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules 
     56   PUBLIC   eos_bn2    ! not use outside eosbn2 module for the moment 
     57   PUBLIC   eos_fzp    ! called by sbcice_... modules 
     58   PUBLIC   eos_pen    ! used for pe diagnostics in trdpen module 
     59   PUBLIC   eos_init   ! called by istate module 
     60 
     61   !                                   !!* Namelist (nameos) * 
     62   INTEGER , PUBLIC ::   nn_eos   = 0   !: = -1/0/1 type of equation of state 
     63 
     64   !                                   !!!  simplified eos coefficients 
     65   REAL(wp) ::   rn_alph0  = 1.67e-4    ! thermal expansion coeff.  
     66   REAL(wp) ::   rn_beta0  = 7.8e-4     ! saline  expansion coeff.  
     67   REAL(wp) ::   rn_gamm0  = 4.4e-6     ! depth expansion coeff.   
     68   REAL(wp) ::   rn_cabbel = 3.0e-2     ! cabbeling coeff.          
     69   REAL(wp) ::   rn_thermo = 1.1e-4     ! thermobaric coeff.       
     70   REAL(wp) ::   offset    = 0          ! offset so that rho(10,35,0)=1027kg/m3 
    7171 
    7272   !! * Substitutions 
     
    7474#  include "vectopt_loop_substitute.h90" 
    7575   !!---------------------------------------------------------------------- 
    76    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     76   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
    7777   !! $Id$ 
    7878   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    8080CONTAINS 
    8181 
    82    SUBROUTINE eos_insitu( pts, prd ) 
    83       !!---------------------------------------------------------------------- 
    84       !!                   ***  ROUTINE eos_insitu  *** 
    85       !! 
    86       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    87       !!       potential temperature and salinity using an equation of state 
    88       !!       defined through the namelist parameter nn_eos. 
    89       !! 
    90       !! ** Method  :   3 cases: 
    91       !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
    92       !!         Check value: rho = 1041.83267 kg/m**3 for p=3000 dbar, 
    93       !!          t = 3 deg celcius, s=35.5 psu 
    94       !!      nn_eos = 0 : standard NEMO equation of state, based on  
    95       !!         Jackett and McDougall (1995) equation of state. 
    96       !!         the in situ density is computed directly as a function of 
    97       !!         potential temperature relative to the surface (the opa t 
    98       !!         variable), salt and pressure (assuming no pressure variation 
    99       !!         along geopotential surfaces, i.e. the pressure p in decibars 
    100       !!         is approximated by the depth in meters. 
    101       !!              prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 
    102       !!         with depth                         z        meters 
    103       !!              potential temperature         t        deg celsius 
    104       !!              salinity                      s        psu 
    105       !!              reference volumic mass        rau0     kg/m**3 
    106       !!              in situ volumic mass          rho      kg/m**3 
    107       !!              in situ density anomalie      prd      no units 
    108       !!         Check value: rho = 1060.93298 kg/m**3 for z=10000 dbar, 
    109       !!          t = 40 deg celcius, s=40 psu 
    110       !!      nn_eos = 1 : Vallis simplified equation of state (Vallis book, 2006) 
    111       !!              prd(t,s,z) =  ( rho(t,s,p) - rau0 ) / rau0 
    112       !!                         = -alpha*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta*(S-S0) + gamma*z 
     82   SUBROUTINE eos_rhd_pot( pts, pdep, prd, prhop ) 
     83      !!---------------------------------------------------------------------- 
     84      !!                   ***  ROUTINE eos_rhd_pot  *** 
     85      !! 
     86      !! ** Purpose :   Compute the density anomaly (rho/rho0-1) from potential 
     87      !!              temperature and salinity and optionally the density referenced  
     88      !!              to the surface using an equation of state defined through 
     89      !!              nn_eos namelist parameter. 
     90      !! 
     91      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 
     92      !!         with   prd    in situ density anomalie     no units 
     93      !!                t      potential temperature        Celsius 
     94      !!                s      salinity                     psu 
     95      !!                z      depth                        meters 
     96      !!                rho    in situ density              kg/m^3 
     97      !!                rau0   reference density            kg/m^3 
     98      !! 
     99      !!     nn_eos = -1 : Jackett and McDougall (1995) equation of state is used for rho(t,s,z). 
     100      !!         Check value: rho = 1041.83267 kg/m^3 for p=3000 dbar, t=3 Celcius, s=35.5 psu 
     101      !! 
     102      !!     nn_eos =  0 : standard NEMO equation of state, based also on Jackett and McDougall (1995)  
     103      !!         but with slightly different coefficients. 
     104      !!         Check value: rho = 1060.93298 kg/m^3 for z=10000 dbar, t=40 Celcius, s=40 psu 
     105      !! 
     106      !!     nn_eos =  1 : simplified equation of state (inspired Vallis 2006) 
     107      !!              prd(t,s,z) = -alpha*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta*(S-S0) + gamma*z 
    113108      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0 
    114109      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 
    115       !!              Vallis (2006) equation: use default values of coefficients  
    116       !!      Note that no boundary condition problem occurs in this routine 
    117       !!      as pts are defined over the whole domain. 
    118       !! 
    119       !! ** Action  :   compute prd , the in situ density (no units) 
     110      !!              Vallis like equation: use default values of coefficients 
     111      !! 
     112      !! ** Action  :   prd   density anomaly (no units) 
    120113      !! 
    121114      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    122115      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
    123116      !!---------------------------------------------------------------------- 
    124       !! 
    125       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    126       !                                                      ! 2 : salinity               [psu] 
    127       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-] 
    128       !! 
     117      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   )           ::   pts    ! pot. temperature and salinity [Celcius,psu] 
     118      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   )           ::   pdep   ! depth of pts data                       [m] 
     119      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::   prd    ! density anomaly                         [-] 
     120      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out), OPTIONAL ::   prhop  ! potential density (ref. to z=0)     [kg/m3] 
     121      ! 
    129122      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     123      INTEGER  ::   ikm1                 ! local integer 
    130124      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
    131125      REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
     
    133127      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
    134128      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
    135       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    136       !!---------------------------------------------------------------------- 
    137  
    138       ! 
    139       IF( nn_timing == 1 ) CALL timing_start('eos') 
     129      !!---------------------------------------------------------------------- 
     130      ! 
     131      IF( nn_timing == 1 )   CALL timing_start('eos_rhd_pot') 
     132      ! 
     133      ikm1 = MAX( SIZE( pts, 3 ) - 1 , 1 )    ! =1 for 2D calculation, jpkm1 otherwise 
    140134      ! 
    141135      SELECT CASE( nn_eos ) 
     
    143137      CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==! 
    144138         ! 
    145          CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    146          ! 
    147 !CDIR NOVERRCHK 
    148          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    149          ! 
    150          DO jk = 1, jpkm1 
     139         DO jk = 1, ikm1 
    151140            DO jj = 1, jpj 
    152141               DO ji = 1, jpi 
    153                   zt = pts  (ji,jj,jk,jp_tem) 
    154                   zs = pts  (ji,jj,jk,jp_sal) 
    155                   zh = fsdept(ji,jj,jk)        ! depth 
    156                   zsr= zws   (ji,jj,jk)        ! square root salinity 
     142                  zt  = pts (ji,jj,jk,jp_tem) 
     143                  zs  = pts (ji,jj,jk,jp_sal) 
     144                  zh  = pdep(ji,jj,jk)                          ! depth 
     145                  zsr = SQRT( ABS( pts(ji,jj,jk,jp_sal) ) )     ! square root salinity 
    157146                  ! 
    158147                  ! compute volumic mass pure water at atm pressure 
     
    182171                  ! 
    183172                  ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 
    184                   prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
    185                      &             - 1._wp  ) * tmask(ji,jj,jk) 
     173                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) * r1_rau0 - 1._wp  ) * tmask(ji,jj,jk) 
     174                  ! 
     175                  IF( PRESENT(prhop) ) THEN                         ! potential density referenced at the surface 
     176                     prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 
     177                  ENDIF 
    186178               END DO 
    187179            END DO 
    188180         END DO 
    189181         ! 
    190          CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    191          ! 
    192182      CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==! 
    193183         ! 
    194          CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    195          ! 
    196 !CDIR NOVERRCHK 
    197          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    198          ! 
    199          DO jk = 1, jpkm1 
     184         DO jk = 1, ikm1 
    200185            DO jj = 1, jpj 
    201186               DO ji = 1, jpi 
    202187                  zt = pts   (ji,jj,jk,jp_tem) 
    203188                  zs = pts   (ji,jj,jk,jp_sal) 
    204                   zh = fsdept(ji,jj,jk)        ! depth 
    205                   zsr= zws   (ji,jj,jk)        ! square root salinity 
     189                  zh = fsdept(ji,jj,jk)                      ! depth 
     190                  zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) )   ! square root salinity 
    206191                  ! 
    207192                  ! compute volumic mass pure water at atm pressure 
     
    230215                  zc = ( zc3*zsr + zc2 )*zs + zc1 
    231216                  ! 
    232                   ! masked in situ density anomaly 
    233                   prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
    234                      &             - 1._wp  ) * tmask(ji,jj,jk) 
     217                  !                                                 ! density anomaly (masked) 
     218                  prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) * r1_rau0 - 1._wp  ) * tmask(ji,jj,jk) 
     219                  ! 
     220                  IF( PRESENT(prhop) ) THEN                         ! potential density referenced at the surface 
     221                     prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 
     222                  ENDIF 
    235223               END DO 
    236224            END DO 
    237225         END DO 
    238226         ! 
    239          CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    240          ! 
    241       CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
    242          DO jk = 1, jpkm1 
     227      CASE( 1 )                !==  simplified EOS  ==! 
     228         DO jk = 1, ikm1 
    243229            DO jj = 1, jpj 
    244230               DO ji = 1, jpi 
    245                   zt = pts   (ji,jj,jk,jp_tem) - 10._wp 
    246                   zs = pts   (ji,jj,jk,jp_sal) - 35._wp 
    247                   zh = fsdept(ji,jj,jk)        ! depth 
    248                   ! masked in situ density anomaly 
    249                   prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          & 
    250                      &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk) 
     231                  zt = pts (ji,jj,jk,jp_tem) - 10._wp 
     232                  zs = pts (ji,jj,jk,jp_sal) - 35._wp 
     233                  zh = pdep(ji,jj,jk) 
     234                  !                                                 ! density anomaly (masked) 
     235                  prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) * zt   & 
     236                     &                     + rn_beta0 * zs + rn_gamm0 * zh                           ) * tmask(ji,jj,jk) 
     237                  ! 
     238                  IF( PRESENT(prhop) ) THEN                         ! potential density referenced at z=0 
     239                     prhop(ji,jj,jk) = rau0 * ( 1._wp                                                  & 
     240                        &         + offset - rn_alph0 * ( 1._wp + rn_cabbel*zt               ) * zt    & 
     241                        &                  + rn_beta0 * zs                                           ) * tmask(ji,jj,jk) 
     242                  ENDIF 
    251243               END DO 
    252244            END DO 
     
    255247      END SELECT 
    256248      ! 
    257       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk ) 
    258       ! 
    259       IF( nn_timing == 1 ) CALL timing_stop('eos') 
    260       ! 
    261    END SUBROUTINE eos_insitu 
    262  
    263  
    264    SUBROUTINE eos_insitu_pot( pts, prd, prhop ) 
    265       !!---------------------------------------------------------------------- 
    266       !!                  ***  ROUTINE eos_insitu_pot  *** 
    267       !! 
    268       !! ** Purpose :   Compute the in situ density (ratio (rho-rau0)/rau0) and the 
    269       !!      potential volumic mass (Kg/m3) from potential temperature and 
    270       !!      salinity fields using an equation of state defined through the 
    271       !!     namelist parameter nn_eos. 
    272       !! 
    273       !! ** Method  :   Same as for eos_insitu 
    274       !! 
    275       !! ** Action  : - prd  , the in situ density (no units) 
    276       !!              - prhop, the potential volumic mass (Kg/m3) 
    277       !! 
    278       !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
    279       !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    280       !!                Brown and Campana, Mon. Weather Rev., 1978 
    281       !!---------------------------------------------------------------------- 
    282       !! 
    283       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius] 
    284       !                                                                ! 2 : salinity               [psu] 
    285       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-] 
    286       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     249      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos_rhd_pot : ', ovlap=1, kdim=jpk ) 
     250      ! 
     251      IF( nn_timing == 1 ) CALL timing_stop('eos_rhd_pot') 
     252      ! 
     253   END SUBROUTINE eos_rhd_pot 
     254 
     255 
     256   SUBROUTINE eos_rab_bn2( pts, pab, pn2 ) 
     257      !!---------------------------------------------------------------------- 
     258      !!                 ***  ROUTINE eos_ab  *** 
     259      !! 
     260      !! ** Purpose :   Calculates the thermal/haline expansion terms at T-points 
     261      !!              and optionally the Brunt-Vaisala frequency square 
     262      !! 
     263      !! ** Method  : - calculates alpha and beta at T-points 
     264      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state 
     265      !!       * nn_eos =  0 : modified Jackett and McDougall (1995) equation of state 
     266      !!       * nn_eos =  1 : Vallis equation of state (Vallis 2006) 
     267      !!                       alpha and beta ratios are returned as 3-D arrays. 
     268      !!              - calculates N^2 through a call to eos_bn2 (if pn2 argument present) 
     269      !! 
     270      !! ** Action  :   pab   partial derivative of density anomaly with respect to T & S at T-points 
     271      !!                pn2   square of the Brunt-Vaisala frequency at w-point  [1/s2] 
     272      !!---------------------------------------------------------------------- 
     273      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   )           ::   pts   ! pot. temperature & salinity [Celcius,psu] 
     274      REAL(wp), DIMENSION(:,:,:,:), INTENT(  out)           ::   pab   ! alpha and beta              [1/Celcius, 1/psu] 
     275      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out), OPTIONAL ::   pn2   ! Brunt-Vaisala frequency**2  [1/s2] 
    287276      ! 
    288277      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     
    292281      REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
    293282      REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
    294       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    295       !!---------------------------------------------------------------------- 
    296       ! 
    297       IF( nn_timing == 1 ) CALL timing_start('eos-p') 
    298       ! 
    299       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    300       ! 
    301       SELECT CASE ( nn_eos ) 
    302       ! 
    303       CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==! 
    304 !CDIR NOVERRCHK 
    305          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     283      REAL(wp) ::   zm_inv, zdm, zdn     !   -      - 
     284      REAL(wp) ::   zdza, zdzb, zdzc     !   -      - 
     285      !!---------------------------------------------------------------------- 
     286      ! 
     287      IF ( nn_timing == 1 )   CALL timing_start('eos_rab_bn2') 
     288      ! 
     289      SELECT CASE ( nn_eos )     !++  partial derivative of density anomaly with respect to T & S  ++! 
     290      ! 
     291      CASE ( -1 )                !==  Jackett and McDougall (1995) formulation  ==! 
    306292         ! 
    307293         DO jk = 1, jpkm1 
     
    310296                  zt = pts   (ji,jj,jk,jp_tem) 
    311297                  zs = pts   (ji,jj,jk,jp_sal) 
    312                   zh = fsdept(ji,jj,jk)        ! depth 
    313                   zsr= zws   (ji,jj,jk)        ! square root salinity 
     298                  zh = fsdept(ji,jj,jk)                        ! depth 
     299                  zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) )   ! square root salinity 
    314300                  ! 
    315301                  ! compute volumic mass pure water at atm pressure 
     
    333319                  zb  = ( zb3*zsr + zb2 ) *zs + zb1 
    334320                  ! 
    335                   zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    336                   zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    337321                  zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
     322                  zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 
     323                  zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 
    338324                  zc  = ( zc3*zsr + zc2 )*zs + zc1 
    339325                  ! 
    340                   ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 
    341                   prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
    342                      &             - 1._wp  ) * tmask(ji,jj,jk) 
    343               prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 
     326                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
     327                  !                                                     ! alpha 
     328                  zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 
     329                  zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   & 
     330                     &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 
     331                  zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   & 
     332                     &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs +    & 
     333                     &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 
     334                  zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     335                     &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     336                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  & 
     337                     &                  -18.19058e-3_wp)*zt+6.793952e-2_wp) 
     338                  zdm = (zdza*zh+zdzb)*zh+zdzc 
     339                  ! 
     340                  pab(ji,jj,jk,jp_tem) = - ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   & 
     341                     &                   * r1_rau0 * tmask(ji,jj,jk) 
     342                  ! 
     343                  zdza = za2                                            ! beta 
     344                  zdzb = 2.220399e-4_wp * zsr + zb2 
     345                  zdzc = 1.5_wp * zc3 * zsr + zc2 
     346                  zdn  = 9.6628e-4_wp * zs + 1.5_wp * zn3 * zsr + zn2 
     347                  zdm  = ( zdza * zh + zdzb ) * zh + zdzc 
     348                  pab(ji,jj,jk,jp_sal) =   ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   & 
     349                     &                   * r1_rau0 * tmask(ji,jj,jk) 
    344350               END DO 
    345351            END DO 
    346352         END DO 
    347353         ! 
    348       CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==! 
    349 !CDIR NOVERRCHK 
    350          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     354      CASE ( 0 )                 !==  modified Jackett and McDougall (1995) formulation  ==! 
    351355         ! 
    352356         DO jk = 1, jpkm1 
     
    355359                  zt = pts   (ji,jj,jk,jp_tem) 
    356360                  zs = pts   (ji,jj,jk,jp_sal) 
    357                   zh = fsdept(ji,jj,jk)        ! depth 
    358                   zsr= zws   (ji,jj,jk)        ! square root salinity 
     361                  zh = fsdept(ji,jj,jk)                      ! depth 
     362                  zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) )   ! square root salinity 
    359363                  ! 
    360364                  ! compute volumic mass pure water at atm pressure 
     
    383387                  zc = ( zc3*zsr + zc2 )*zs + zc1 
    384388                  ! 
    385                   ! masked in situ density anomaly 
    386                   prd(ji,jj,jk) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
    387                      &             - 1._wp  ) * tmask(ji,jj,jk) 
    388               prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 
    389                END DO 
    390             END DO 
    391          END DO 
    392          ! 
    393       CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
    394          DO jk = 1, jpkm1 
    395             DO jj = 1, jpj 
    396                DO ji = 1, jpi 
    397                   zt = pts   (ji,jj,jk,jp_tem) - 10._wp 
    398                   zs = pts   (ji,jj,jk,jp_sal) - 35._wp 
    399                   zh = fsdept(ji,jj,jk)        ! depth 
    400                   ! masked in situ density anomaly 
    401                   prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          & 
    402                      &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk) 
    403                   ! masked in situ density anomaly 
    404                  prhop(ji,jj,jk) = ( 1.0_wp + offset - rn_alph0 * ( 1._wp + rn_cabbel*zt ) *zt + rn_beta0 *zs )  & 
    405                      &              * rau0 * tmask(ji,jj,jk) 
    406                  ! 
    407                END DO 
    408             END DO 
    409          END DO 
    410          ! 
    411       END SELECT 
    412       ! 
    413       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 
    414       ! 
    415       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    416       ! 
    417       IF( nn_timing == 1 ) CALL timing_stop('eos-p') 
    418       ! 
    419    END SUBROUTINE eos_insitu_pot 
    420  
    421  
    422    SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 
    423       !!---------------------------------------------------------------------- 
    424       !!                  ***  ROUTINE eos_insitu_2d  *** 
    425       !! 
    426       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    427       !!      potential temperature and salinity using an equation of state 
    428       !!      defined through the namelist parameter nn_eos. * 2D field case 
    429       !! 
    430       !! ** Method  :   Same as for eos_insitu 
    431       !! 
    432       !! ** Action  : - prd , the in situ density (no units) 
    433       !! 
    434       !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
    435       !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    436       !!---------------------------------------------------------------------- 
    437       !! 
    438       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    439       !                                                           ! 2 : salinity               [psu] 
    440       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m] 
    441       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
    442       !! 
    443       INTEGER  ::   ji, jj              ! dummy loop indices 
    444       REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
    445       REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
    446       REAL(wp) ::   zn,  za1, za2, za    !   -      - 
    447       REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
    448       REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
    449       REAL(wp), POINTER, DIMENSION(:,:) :: zws 
    450       !!---------------------------------------------------------------------- 
    451       ! 
    452       IF( nn_timing == 1 ) CALL timing_start('eos2d') 
    453       ! 
    454       CALL wrk_alloc( jpi, jpj, zws ) 
    455       ! 
    456  
    457       prd(:,:) = 0._wp 
    458  
    459       SELECT CASE( nn_eos ) 
    460       ! 
    461       CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==! 
    462       ! 
    463 !CDIR NOVERRCHK 
    464          DO jj = 1, jpjm1 
    465 !CDIR NOVERRCHK 
    466             DO ji = 1, fs_jpim1   ! vector opt. 
    467                zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 
    468             END DO 
    469          END DO 
    470          DO jj = 1, jpjm1 
    471             DO ji = 1, fs_jpim1   ! vector opt. 
    472                zt    = pts  (ji,jj,jp_tem)     ! interpolated T 
    473                zs    = pts  (ji,jj,jp_sal)     ! interpolated S 
    474                zsr   = zws  (ji,jj)            ! square root of interpolated S 
    475                zh    = pdep (ji,jj)            ! depth at the partial step level 
    476                ! 
    477                ! compute volumic mass pure water at atm pressure 
    478                zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    479                   &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    480                ! seawater volumic mass atm pressure 
    481                zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    482                   &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    483                zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    484                zn4= 4.8314e-4_wp 
    485                zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    486                ! 
    487                ! add the compression terms 
    488                za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    489                za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    490                za = za1 + za2 * zs 
    491                ! 
    492                zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
    493                zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
    494                zb3 = 1.480266e-4_wp 
    495                zb  = ( zb3*zsr + zb2 ) *zs + zb1 
    496                ! 
    497                zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
    498                zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
    499                zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
    500                zc  = ( zc3*zsr + zc2 )*zs + zc1 
    501                ! 
    502                ! masked in situ density anomaly 
    503                prd(ji,jj) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
    504                   &             - 1._wp  ) * tmask(ji,jj,1) 
    505             END DO 
    506          END DO 
    507          ! 
    508       CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==! 
    509       ! 
    510 !CDIR NOVERRCHK 
    511          DO jj = 1, jpjm1 
    512 !CDIR NOVERRCHK 
    513             DO ji = 1, fs_jpim1   ! vector opt. 
    514                zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 
    515             END DO 
    516          END DO 
    517          DO jj = 1, jpjm1 
    518             DO ji = 1, fs_jpim1   ! vector opt. 
    519                zt    = pts  (ji,jj,jp_tem)     ! interpolated T 
    520                zs    = pts  (ji,jj,jp_sal)     ! interpolated S 
    521                zsr   = zws  (ji,jj)            ! square root of interpolated S 
    522                zh    = pdep (ji,jj)            ! depth at the partial step level 
    523                ! 
    524                ! compute volumic mass pure water at atm pressure 
    525                zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    526                   &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    527                ! seawater volumic mass atm pressure 
    528                zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    529                   &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    530                zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    531                zn4= 4.8314e-4_wp 
    532                zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    533                ! 
    534                ! add the compression terms 
    535                za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    536                za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    537                za = za1 + za2 * zs 
    538                ! 
    539                zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
    540                zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
    541                zb3= 2.042967e-2_wp 
    542                zb = ( zb3*zsr + zb2 ) *zs + zb1 
    543                ! 
    544                zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    545                zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
    546                zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
    547                zc = ( zc3*zsr + zc2 )*zs + zc1 
    548                ! 
    549                ! masked in situ density anomaly 
    550                prd(ji,jj) = (  zn * (  1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) )  / rau0   & 
    551                   &             - 1._wp  ) * tmask(ji,jj,1) 
    552             END DO 
    553          END DO 
    554          ! 
    555       CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
    556          DO jj = 1, jpjm1 
    557             DO ji = 1, fs_jpim1   ! vector opt. 
    558                zt    = pts  (ji,jj,jp_tem) - 10._wp   ! interpolated T 
    559                zs    = pts  (ji,jj,jp_sal) - 35._wp   ! interpolated S 
    560                zh    = pdep (ji,jj)            ! depth at the partial step level 
    561                ! masked in situ density anomaly 
    562                prd(ji,jj) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt  &                                           & 
    563                      &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,1) 
    564                ! 
    565             END DO 
    566          END DO 
    567          ! 
    568       END SELECT 
    569  
    570       IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 
    571       ! 
    572       CALL wrk_dealloc( jpi, jpj, zws ) 
    573       ! 
    574       IF( nn_timing == 1 ) CALL timing_stop('eos2d') 
    575       ! 
    576    END SUBROUTINE eos_insitu_2d 
    577  
    578  
    579  
    580    SUBROUTINE eos_ab( pts, pab ) 
    581       !!---------------------------------------------------------------------- 
    582       !!                 ***  ROUTINE eos_ab  *** 
    583       !! 
    584       !! ** Purpose :   Calculates the in situ thermal/haline expansion terms at T-points 
    585       !! 
    586       !! ** Method  :   calculates alpha and beta at T-points 
    587       !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state 
    588       !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state 
    589       !!       * nn_eos = 1  : Vallis equation of state (Vallis 2006) 
    590       !!                       alpha and beta ratios are returned as 3-D arrays. 
    591       !! 
    592       !! ** Action  : - pab : partial derivative of in situ density with respect to T & S at T-points 
    593       !!---------------------------------------------------------------------- 
    594       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity 
    595       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! alpha and beta 
    596       ! 
    597       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    598       REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
    599       REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
    600       REAL(wp) ::   zn,  za1, za2, za    !   -      - 
    601       REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
    602       REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
    603       REAL(wp) ::   zm_inv, zdm, zdn     !   -      - 
    604       REAL(wp) ::   zdza, zdzb, zdzc     !   -      - 
    605       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    606       !!---------------------------------------------------------------------- 
    607       ! 
    608       IF ( nn_timing == 1 )   CALL timing_start('eos_ab') 
    609       ! 
    610       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    611       ! 
    612       SELECT CASE ( nn_eos ) 
    613       ! 
    614       CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
    615          ! 
    616          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    617          DO jk = 1, jpkm1 
    618             DO jj = 1, jpj 
    619                DO ji = 1, jpi 
    620                   zt = pts   (ji,jj,jk,jp_tem) 
    621                   zs = pts   (ji,jj,jk,jp_sal) 
    622                   zh = fsdept(ji,jj,jk)        ! depth 
    623                   zsr= zws   (ji,jj,jk)        ! square root salinity 
    624                   ! 
    625                   ! compute volumic mass pure water at atm pressure 
    626                   zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    627                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    628                   ! seawater volumic mass atm pressure 
    629                   zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    630                      &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    631                   zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    632                   zn4= 4.8314e-4_wp 
    633                   zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    634                   ! 
    635                   ! add the compression terms 
    636                   za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    637                   za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    638                   za = za1 + za2 * zs 
    639                   ! 
    640                   zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
    641                   zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
    642                   zb3 = 1.480266e-4_wp 
    643                   zb  = ( zb3*zsr + zb2 ) *zs + zb1 
    644                   ! 
    645                   zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
    646                   zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 
    647                   zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 
    648                   zc  = ( zc3*zsr + zc2 )*zs + zc1 
    649                   ! 
    650389                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
    651               ! alpha 
    652                   zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 
    653                   zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   & 
    654                      &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 
    655                   zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   & 
    656                      &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + & 
    657                      &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 
    658                       
    659                   zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    660                      &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    661                      &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  & 
    662                      &                  -18.19058e-3_wp)*zt+6.793952e-2_wp) 
    663                   zdm = (zdza*zh+zdzb)*zh+zdzc 
    664                   ! 
    665                   pab(ji,jj,jk,jp_tem) =  - ( & 
    666                      &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
    667                      &   / rau0 * tmask(ji,jj,jk) 
    668                   ! 
    669                   ! beta 
    670                   zdza  = za2 
    671                   zdzb  = 2.220399e-4_wp*zsr+zb2 
    672                   zdzc  = 1.5_wp*zc3*zsr+zc2 
    673                   zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
    674                   zdm   = (zdza*zh+zdzb)*zh+zdzc 
    675                   pab(ji,jj,jk,jp_sal) =    ( & 
    676                      &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
    677                      &   / rau0 * tmask(ji,jj,jk) 
    678                   ! 
    679                END DO 
    680             END DO 
    681          END DO 
    682          ! 
    683       CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
    684          ! 
    685          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    686          DO jk = 1, jpkm1 
    687             DO jj = 1, jpj 
    688                DO ji = 1, jpi 
    689                   zt = pts   (ji,jj,jk,jp_tem) 
    690                   zs = pts   (ji,jj,jk,jp_sal) 
    691                   zh = fsdept(ji,jj,jk)        ! depth 
    692                   zsr= zws   (ji,jj,jk)        ! square root salinity 
    693                   ! 
    694                   ! compute volumic mass pure water at atm pressure 
    695                   zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    696                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    697                   ! seawater volumic mass atm pressure 
    698                   zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    699                      &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    700                   zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    701                   zn4= 4.8314e-4_wp 
    702                   zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    703                   ! 
    704                   ! add the compression terms 
    705                   za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    706                   za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    707                   za = za1 + za2 * zs 
    708                   ! 
    709                   zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
    710                   zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
    711                   zb3= 2.042967e-2_wp 
    712                   zb = ( zb3*zsr + zb2 ) *zs + zb1 
    713                   ! 
    714                   zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    715                   zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
    716                   zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
    717                   zc = ( zc3*zsr + zc2 )*zs + zc1 
    718                   ! 
    719                   zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
    720               ! alpha 
     390                  !                                                     ! alpha 
    721391                  zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    722392                  zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp 
     
    727397                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    728398                  zdm = (zdza*zh+zdzb)*zh+zdzc 
    729                   ! 
    730                   pab(ji,jj,jk,jp_tem) =  - ( & 
    731                      &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
    732                      &   / rau0 * tmask(ji,jj,jk) 
    733                   ! beta 
    734                   zdza  = za2 
     399                  pab(ji,jj,jk,jp_tem) = - ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   & 
     400                     &                   * r1_rau0 * tmask(ji,jj,jk) 
     401                  ! 
     402                  zdza  = za2                                           ! beta 
    735403                  zdzb  = 3.0644505e-2_wp*zsr+zb2 
    736404                  zdzc  = 1.5_wp*zc3*zsr+zc2 
    737405                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
    738406                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
    739                   pab(ji,jj,jk,jp_sal) =    ( & 
    740                      &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
    741                      &   / rau0 * tmask(ji,jj,jk) 
    742                   ! 
     407                  pab(ji,jj,jk,jp_sal) =   ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   & 
     408                     &                   * r1_rau0 * tmask(ji,jj,jk) 
    743409               END DO 
    744410            END DO 
    745411         END DO 
    746412         ! 
    747       CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
     413      CASE( 1 )                  !== simplified EOS  ==! 
    748414         DO jk = 1, jpkm1 
    749415            DO jj = 1, jpj 
    750416               DO ji = 1, jpi 
    751                   zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature anomaly (t-T0) 
    752                   zh = fsdept(ji,jj,jk)               ! depth in meters at t-point 
    753                   ! 
     417                  zt = pts (ji,jj,jk,jp_tem) - 10._wp  ! pot. temperature anomaly (t-T0) 
     418                  zh = fsdept(ji,jj,jk)                ! depth in meters at t-point 
     419                  !                                    ! alpha 
    754420                  pab(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*rn_cabbel*zt + rn_thermo*zh ) * tmask(ji,jj,jk)   ! alpha 
     421                  !                                    ! beta = constant: set one for all during the initialisation 
    755422               END DO 
    756423            END DO 
    757424         END DO 
    758          pab(:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! beta 
    759          ! 
    760       CASE DEFAULT 
    761          IF(lwp) WRITE(numout,cform_err) 
    762          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
    763          nstop = nstop + 1 
    764425         ! 
    765426      END SELECT 
    766427      ! 
    767       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    768       ! 
    769       IF( nn_timing == 1 )   CALL timing_stop('eos_ab') 
    770       ! 
    771    END SUBROUTINE eos_ab 
    772  
    773  
    774  
    775    SUBROUTINE eos_insitu_ab( pts, prd, pab ) 
    776       !!---------------------------------------------------------------------- 
    777       !!                 ***  ROUTINE eos_insitu_ab  *** 
    778       !! 
    779       !! ** Purpose :   Calculates the in situ thermal/haline expansion terms 
    780       !!                  and the insitu density anomaly at T-points 
    781       !! 
    782       !! ** Method  :   calculates rhd, alpha, beta at T-points 
    783       !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state 
    784       !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state 
    785       !!       * nn_eos = 1  : Vallis equation of state (Vallis 2006) 
    786       !!                       alpha and beta ratios are returned as 3-D arrays. 
    787       !! 
    788       !! ** Action  : - pab : thermal and haline expansion ratios at T-points 
    789       !!              - prd , the density anomaly (no units) 
     428      !                          !++  Brunt-Vaisala frequency  (only for 3D computation)  ++! 
     429      IF( PRESENT( pn2 ) )   CALL eos_bn2( pts, pab, pn2 )       ! N^2 
     430      ! 
     431      IF( nn_timing == 1 )   CALL timing_stop('eos_rab_bn2') 
     432      ! 
     433   END SUBROUTINE eos_rab_bn2 
     434 
     435 
     436   SUBROUTINE eos_bn2( pts, pab, pn2 ) 
     437      !!---------------------------------------------------------------------- 
     438      !!                  ***  ROUTINE eos_bn2  *** 
     439      !! 
     440      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the  
     441      !!                time-step of the input arguments 
     442      !! 
     443      !! ** Method  :   pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w 
     444      !!      where alpha and beta have been vertically interpolated from T- to w-point. 
     445      !!      N.B. N^2 is set one for all to zero at jk=1 in istate module. 
     446      !! 
     447      !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point  
    790448      !! 
    791449      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
    792450      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    793       !!---------------------------------------------------------------------- 
    794       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity    [Celcius,psu] 
    795       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   prd   ! density anomaly                [-] 
    796       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! partial derivative of density anomaly 
    797       !                                                               ! with respect to T and S        [1/Celcius,1/psu] 
    798       ! 
    799       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    800       REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
    801       REAL(wp) ::   zn1, zn2, zn3, zn4   !   -      - 
    802       REAL(wp) ::   zn,  za1, za2, za    !   -      - 
    803       REAL(wp) ::   zb1, zb2, zb3, zb    !   -      - 
    804       REAL(wp) ::   zc1, zc2, zc3, zc    !   -      - 
    805       REAL(wp) ::   zm_inv, zdm, zdn     !   -      - 
    806       REAL(wp) ::   zdza, zdzb, zdzc     !   -      - 
    807       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    808       !!---------------------------------------------------------------------- 
    809       ! 
    810       IF ( nn_timing == 1 )   CALL timing_start('eos_insitu_ab') 
    811       ! 
    812       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    813       ! 
    814       SELECT CASE ( nn_eos ) 
    815       ! 
    816       CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
    817          ! 
    818          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    819          DO jk = 1, jpkm1 
    820             DO jj = 1, jpj 
    821                DO ji = 1, jpi 
    822                   zt = pts   (ji,jj,jk,jp_tem) 
    823                   zs = pts   (ji,jj,jk,jp_sal) 
    824                   zh = fsdept(ji,jj,jk)        ! depth 
    825                   zsr= zws   (ji,jj,jk)        ! square root salinity 
    826                   ! 
    827                   ! compute volumic mass pure water at atm pressure 
    828                   zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    829                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    830                   ! seawater volumic mass atm pressure 
    831                   zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    832                      &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    833                   zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    834                   zn4= 4.8314e-4_wp 
    835                   zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    836                   ! 
    837                   ! add the compression terms 
    838                   za1= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
    839                   za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
    840                   za = za1 + za2 * zs 
    841                   ! 
    842                   zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 
    843                   zb2 =   ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 
    844                   zb3 = 1.480266e-4_wp 
    845                   zb  = ( zb3*zsr + zb2 ) *zs + zb1 
    846                   ! 
    847                   zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 
    848                   zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 
    849                   zc3 =   (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 
    850                   zc  = ( zc3*zsr + zc2 )*zs + zc1 
    851                   ! 
    852                   zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
    853                   ! 
    854                   ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 
    855                   prd(ji,jj,jk) = (  zn * (  1.0_wp + zh * zm_inv ) / rau0  - 1._wp  ) * tmask(ji,jj,jk) 
    856                   ! 
    857               ! alpha 
    858                   zdza  = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 
    859                   zdzb  = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs +   & 
    860                      &        ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 
    861                   zdzc  = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr +   & 
    862                      &       (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + & 
    863                      &       ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 
    864                       
    865                   zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    866                      &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    867                      &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt  & 
    868                      &                  -18.19058e-3_wp)*zt+6.793952e-2_wp) 
    869                   zdm = (zdza*zh+zdzb)*zh+zdzc 
    870                   ! 
    871                   pab(ji,jj,jk,jp_tem) =  - ( & 
    872                      &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
    873                      &   / rau0 * tmask(ji,jj,jk) 
    874                   ! 
    875                   ! beta 
    876                   zdza  = za2 
    877                   zdzb  = 2.220399e-4_wp*zsr+zb2 
    878                   zdzc  = 1.5_wp*zc3*zsr+zc2 
    879                   zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
    880                   zdm   = (zdza*zh+zdzb)*zh+zdzc 
    881                   pab(ji,jj,jk,jp_sal) =    ( & 
    882                      &   zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 
    883                      &   / rau0 * tmask(ji,jj,jk) 
    884                   ! 
    885                END DO 
    886             END DO 
    887          END DO 
    888          ! 
    889       CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
    890          ! 
    891          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    892          DO jk = 1, jpkm1 
    893             DO jj = 1, jpj 
    894                DO ji = 1, jpi 
    895                   zt = pts   (ji,jj,jk,jp_tem) 
    896                   zs = pts   (ji,jj,jk,jp_sal) 
    897                   zh = fsdept(ji,jj,jk)        ! depth 
    898                   zsr= zws   (ji,jj,jk)        ! square root salinity 
    899                   ! 
    900                   ! compute volumic mass pure water at atm pressure 
    901                   zn1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
    902                      &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
    903                   ! seawater volumic mass atm pressure 
    904                   zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
    905                      &                      -4.0899e-3_wp ) *zt+0.824493_wp 
    906                   zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
    907                   zn4= 4.8314e-4_wp 
    908                   zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 
    909                   ! 
    910                   ! add the compression terms 
    911                   za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    912                   za1= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    913                   za = za1 + za2 * zs 
    914                   ! 
    915                   zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 
    916                   zb2=   (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 
    917                   zb3= 2.042967e-2_wp 
    918                   zb = ( zb3*zsr + zb2 ) *zs + zb1 
    919                   ! 
    920                   zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
    921                   zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
    922                   zc3=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
    923                   zc = ( zc3*zsr + zc2 )*zs + zc1 
    924                   ! 
    925                   zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
    926                   ! 
    927                   ! masked in situ density anomaly (rho/rho0 - 1). Important: rau0=1035, should be 1027! 
    928                   prd(ji,jj,jk) = (  zn * (  1.0_wp + zh * zm_inv ) * r1_rau0  - 1._wp  ) * tmask(ji,jj,jk) 
    929                   ! 
    930                   !                                 ! alpha 
    931                   zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    932                   zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp 
    933                   zdzc  = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
    934                      &    +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
    935                   zdn   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
    936                      &    +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
    937                      &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    938                   zdm = (zdza*zh+zdzb)*zh+zdzc 
    939                   ! 
    940                   pab(ji,jj,jk,jp_tem) =  - (  zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   & 
    941                      &                    * r1_rau0 * tmask(ji,jj,jk) 
    942                   ! 
    943                   !                                 ! beta 
    944                   zdza  = za2 
    945                   zdzb  = 3.0644505e-2_wp*zsr+zb2 
    946                   zdzc  = 1.5_wp*zc3*zsr+zc2 
    947                   zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
    948                   zdm   = (zdza*zh+zdzb)*zh+zdzc 
    949                   pab(ji,jj,jk,jp_sal) = (  zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv )   & 
    950                      &                 * r1_rau0 * tmask(ji,jj,jk) 
    951                END DO 
    952             END DO 
    953          END DO 
    954          ! 
    955       CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
    956          DO jk = 1, jpkm1 
    957             DO jj = 1, jpj 
    958                DO ji = 1, jpi 
    959                   zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-T0) 
    960                   zh = fsdept(ji,jj,jk)                                                ! depth in meters  at t-point 
    961                   ! masked in situ density anomaly 
    962                   prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt   &                                          & 
    963                      &   + rn_beta0 * zs + rn_gamm0 * zh  ) * tmask(ji,jj,jk) 
    964                   ! alpha 
    965                   pab(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*rn_cabbel*zt + rn_thermo*zh )  & 
    966                       &                       * tmask(ji,jj,jk)   ! alpha 
    967                   ! 
    968                END DO 
    969             END DO 
    970          END DO 
    971          ! beta 
    972          pab (:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! beta 
    973          ! 
    974       CASE DEFAULT 
    975          IF(lwp) WRITE(numout,cform_err) 
    976          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
    977          nstop = nstop + 1 
    978          ! 
    979       END SELECT 
    980       ! 
    981       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    982       ! 
    983       IF( nn_timing == 1 ) CALL timing_stop('eos_insitu_ab') 
    984       ! 
    985    END SUBROUTINE eos_insitu_ab 
    986  
     451      !!                McDougall, J. Phys. Oceano., 1987 
     452      !!---------------------------------------------------------------------- 
     453      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu] 
     454      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::  pab   ! partial derivative of in situ  
     455      !                                                              ! density with respect to T and S [1/Celcius,1/psu] 
     456      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency    [1/s^2] 
     457      !! 
     458      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     459      REAL(wp) ::   zaw, zbw, zrw   ! local scalars 
     460      !!---------------------------------------------------------------------- 
     461      ! 
     462      IF( nn_timing == 1 ) CALL timing_start('eos_bn2') 
     463      ! 
     464      ! pn2 : interior points only (2=< jk =< jpkm1 ) 
     465      ! -------------------------- 
     466      ! 
     467      DO jk = 2, jpkm1 
     468         DO jj = 1, jpj 
     469            DO ji = 1, jpi 
     470               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   & 
     471                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )  
     472               ! 
     473               zaw = (  pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  ) * tmask(ji,jj,jk) 
     474               zbw = (  pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw  ) * tmask(ji,jj,jk) 
     475               ! 
     476               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
     477                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  ) / fse3w(ji,jj,jk) 
     478            END DO 
     479         END DO 
     480      END DO 
     481      ! 
     482      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk ) 
     483      ! 
     484      IF( nn_timing == 1 )   CALL timing_stop('eos_bn2') 
     485      ! 
     486   END SUBROUTINE eos_bn2 
     487 
     488 
     489   FUNCTION eos_fzp( psal, pdep ) RESULT( ptf ) 
     490      !!---------------------------------------------------------------------- 
     491      !!                 ***  ROUTINE eos_fzp  *** 
     492      !! 
     493      !! ** Purpose :   Compute the freezing point temperature [Celcius] 
     494      !! 
     495      !! ** Method  :   UNESCO freezing point (ptf) in Celcius is given by 
     496      !!       ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z 
     497      !!       checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m 
     498      !! 
     499      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978 
     500      !!---------------------------------------------------------------------- 
     501      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu] 
     502      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m] 
     503      ! Leave result array automatic rather than making explicitly allocated 
     504      REAL(wp), DIMENSION(jpi,jpj) ::   ptf   ! freezing temperature [Celcius] 
     505      !!---------------------------------------------------------------------- 
     506      ! 
     507      ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
     508         &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:) 
     509      ! 
     510      IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 
     511      ! 
     512   END FUNCTION eos_fzp 
    987513 
    988514 
     
    1027553      REAL(wp) ::   zddelta, zdAp        !   -      - 
    1028554      REAL(wp) ::   zdlogBp, zdCp, zdP   !   -      - 
    1029       REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    1030555      !!---------------------------------------------------------------------- 
    1031556      ! 
    1032557      IF ( nn_timing == 1 )   CALL timing_start('eos_pen') 
    1033558      ! 
    1034       CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    1035       ! 
    1036559      SELECT CASE ( nn_eos ) 
    1037560      ! 
    1038561      CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
    1039562         ! 
    1040          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    1041563         DO jk = 1, jpkm1 
    1042564            DO jj = 1, jpj 
     
    1044566                  zt = pts   (ji,jj,jk,jp_tem) 
    1045567                  zs = pts   (ji,jj,jk,jp_sal) 
    1046                   zh = fsdept(ji,jj,jk)        ! depth (Important: zh must be different from 0) 
    1047                   zsr= zws   (ji,jj,jk)        ! square root salinity 
     568                  zh = fsdept(ji,jj,jk)                      ! depth (Important: zh must be different from 0) 
     569                  zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) )   ! square root salinity 
    1048570                  ! 
    1049571                  ! compute volumic mass pure water at atm pressure 
     
    1100622                  ! 
    1101623                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
    1102                   zsgn   = SIGN( 1. , zdelta ) 
    1103                   zsdelta = SQRT( zsgn * zdelta ) 
     624                  zsdelta = SQRT(  ABS( zdelta )  ) 
    1104625                  zAp     = - zb / za / zsdelta 
    1105626                  zlogBp  = - LOG( zc * zm_inv ) 
    1106627                  zCp     = zh * zsdelta / (2._wp*zc+zh*zb) 
    1107                   IF( zsgn > 0. ) THEN ;   zatCp = ATAN(zCp); 
    1108                   ELSE  ;    zatCp = ATANH(zCp) 
    1109                   ENDIF 
     628                  zsgn    = SIGN( 1._wp , zdelta ) 
     629                  zatCp   = 0.5_wp * (  ATAN( zCp ) * (1._wp + zsgn) + ATANH( zCp ) * (1._wp - zsgn)  ) 
    1110630                  zP      = ( 0.5_wp*zlogBp/za + zAp * zatCp ) / zh 
    1111631                  ! 
    1112                   ! compute potential energy 
    1113                   pen(ji,jj,jk)  = ( zn - rau0 + zn * zP ) / rau0 * tmask(ji,jj,jk) 
    1114                   ! 
    1115               ! alphaPE 
     632                  !                              ! potential energy 
     633                  pen(ji,jj,jk)  = ( zn - rau0 + zn * zP ) * r1_rau0 * tmask(ji,jj,jk) 
     634                  ! 
     635                  !                              ! alphaPE 
    1116636                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 
    1117637                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta) 
     
    1121641                     &       + zAp*zdCp / ( 1._wp + zsgn*zCp*zCp )  ) / zh 
    1122642                  ! 
    1123                   pab_pe(ji,jj,jk,jp_tem)    = - ( zdn * (1._wp + zP ) + zn * zdP ) / rau0 * tmask(ji,jj,jk) 
    1124                   ! dds 
    1125                   zdza  = za2 
     643                  pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) * r1_rau0 * tmask(ji,jj,jk) 
     644                  ! 
     645                  zdza  = za2                    ! dds 
    1126646                  zdzb  = 2.220399e-4_wp*zsr+zb2 
    1127647                  zdzc  = 1.5_wp*zc3*zsr+zc2 
    1128648                  zdn   = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 
    1129649                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
    1130                   ! betaPE 
     650                  !                              ! betaPE 
    1131651                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 
    1132652                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta) 
     
    1136656                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 
    1137657                  ! 
    1138                   pab_pe(ji,jj,jk,jp_sal)    = ( zdn * (1._wp + zP ) + zn * zdP )    & 
    1139                                                     &           / rau0 * tmask(ji,jj,jk) 
    1140                   ! 
     658                  pab_pe(ji,jj,jk,jp_sal) =   ( zdn * (1._wp + zP ) + zn * zdP ) * r1_rau0 * tmask(ji,jj,jk) 
    1141659               END DO 
    1142660            END DO 
     
    1145663      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
    1146664         ! 
    1147          zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    1148665         DO jk = 1, jpkm1 
    1149666            DO jj = 1, jpj 
    1150667               DO ji = 1, jpi 
    1151                   ! 
    1152668                  zt = pts   (ji,jj,jk,jp_tem) 
    1153669                  zs = pts   (ji,jj,jk,jp_sal) 
    1154                   zh = fsdept(ji,jj,jk)        ! depth (Important: zh must be different from 0) 
    1155                   zsr= zws   (ji,jj,jk)        ! square root salinity 
     670                  zh = fsdept(ji,jj,jk)                      ! depth (Important: zh must be different from 0) 
     671                  zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) )   ! square root salinity 
    1156672                  ! 
    1157673                  ! compute volumic mass pure water at atm pressure 
     
    1182698                  zdelta = 4._wp * za * zc - zb * zb 
    1183699                  ! 
    1184                   ! dds 
    1185                   zdza  = za2 
     700                  zdza  = za2                    ! dds 
    1186701                  zdzb  = 3.0644505e-2_wp*zsr+zb2 
    1187702                  zdzc  = 1.5_wp*zc3*zsr+zc2 
     
    1189704                  zdm   = (zdza*zh+zdzb)*zh+zdzc 
    1190705                  ! stability condition for the calculation (add a small perturbation on S if needed) 
    1191                   IF( ABS(zdelta) < 1.e-5 .OR. ABS(za)<1.e-9 ) THEN ; !stability condition for the calculation 
     706                  IF( ABS(zdelta) < 1.e-5 .OR. ABS(za) < 1.e-9 ) THEN  !stability condition for the calculation 
    1192707                     zs = zs + 1.e-3_wp 
    1193708                     zsr= zsr + 0.5e-3_wp/zsr 
     
    1199714                     ! 
    1200715                     zdelta = 4._wp * za * zc - zb * zb 
    1201                      ! 
    1202716                  ENDIF 
    1203717                  ! 
    1204718                  zm_inv  = 1._wp / ( ( za*zh + zb )*zh + zc ) 
    1205                   zsgn   = SIGN( 1. , zdelta ) 
    1206                   zsdelta = SQRT( zsgn * zdelta ) 
     719                  zsdelta = SQRT(  ABS( zdelta )  ) 
    1207720                  zAp     = - zb / za / zsdelta 
    1208721                  zlogBp  = - LOG( zc * zm_inv ) 
    1209722                  zCp     = zh * zsdelta / (2._wp*zc+zh*zb) 
    1210                   IF( zsgn > 0. ) THEN ;   zatCp = ATAN(zCp); 
    1211                   ELSE  ;    zatCp = ATANH(zCp) 
    1212                   ENDIF 
     723                  zsgn    = SIGN( 1._wp , zdelta ) 
     724                  zatCp   = 0.5_wp * (  ATAN( zCp ) * (1._wp + zsgn) + ATANH( zCp ) * (1._wp - zsgn)  ) 
    1213725                  zP      = ( 0.5_wp*zlogBp/za + zAp * zatCp ) / zh 
    1214726                  ! 
    1215                   ! compute potential energy 
    1216                   pen(ji,jj,jk)  = ( zn - rau0 + zn * zP ) / rau0 * tmask(ji,jj,jk) 
    1217                   ! 
    1218                   ! betaPE 
     727                  !                               ! potential energy 
     728                  pen(ji,jj,jk)  = ( zn - rau0 + zn * zP ) * r1_rau0 * tmask(ji,jj,jk) 
     729                  ! 
     730                  !                               ! betaPE 
    1219731                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 
    1220732                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta) 
     
    1224736                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 
    1225737                  ! 
    1226                   pab_pe(ji,jj,jk,jp_sal)    = ( zdn * (1._wp + zP ) + zn * zdP )    & 
    1227                                                     &           / rau0 * tmask(ji,jj,jk) 
    1228                   ! 
    1229               ! ddt 
     738                  pab_pe(ji,jj,jk,jp_sal) = ( zdn * (1._wp + zP ) + zn * zdP ) * r1_rau0 * tmask(ji,jj,jk) 
     739                  ! 
     740                  !                               ! ddt 
    1230741                  zdza  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
    1231742                  zdzb  = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp 
     
    1236747                     &    +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
    1237748                  zdm = (zdza*zh+zdzb)*zh+zdzc 
    1238               ! alphaPE 
     749                  !                               ! alphaPE 
    1239750                  zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 
    1240751                  zdAp    = zAp * (zdzb/zb - zdza/za - zddelta) 
     
    1244755                     &          zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 
    1245756                  ! 
    1246                   pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) / rau0 * tmask(ji,jj,jk) 
    1247                   ! 
     757                  pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) * r1_rau0 * tmask(ji,jj,jk) 
    1248758               END DO 
    1249759            END DO 
     
    1254764            DO jj = 1, jpj 
    1255765               DO ji = 1, jpi 
    1256                   zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-T0) 
    1257                   zh = fsdept(ji,jj,jk)                                                ! depth in meters  at t-point 
    1258                   pab_pe(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*(rn_cabbel*zt + rn_thermo*zh) )  & 
    1259                       &          * tmask(ji,jj,jk)                                     ! alphaPE 
     766                  zt = pts(ji,jj,jk,jp_tem) - 10._wp   ! potential temperature (t-T0) 
     767                  zh = fsdept(ji,jj,jk)                ! depth in meters  at t-point 
     768                  !                                    ! alphaPE 
     769                  pab_pe(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*(rn_cabbel*zt + rn_thermo*zh) ) * tmask(ji,jj,jk) 
     770                  !                                    ! Potential Energy 
    1260771                  pen(ji,jj,jk) = ( - rn_alph0 * ( 1._wp + rn_cabbel*zt + 0.5_wp*rn_thermo*zh ) *zt   &                                          & 
    1261772                      &          + rn_beta0 * zs + 0.5_wp*rn_gamm0 * zh  ) * tmask(ji,jj,jk) 
    1262                   ! 
    1263773               END DO 
    1264774            END DO 
     
    1273783      END SELECT 
    1274784      ! 
    1275       CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
    1276       ! 
    1277       IF( nn_timing == 1 ) CALL timing_stop('eos_pen') 
     785      IF( nn_timing == 1 )   CALL timing_stop('eos_pen') 
    1278786      ! 
    1279787   END SUBROUTINE eos_pen 
    1280  
    1281  
    1282  
    1283    SUBROUTINE eos_bn2_ab( pts, pn2 ) 
    1284       !!---------------------------------------------------------------------- 
    1285       !!                  ***  ROUTINE eos_bn2_ab  *** 
    1286       !! 
    1287       !! ** Purpose :   Update alpha/beta then compute the local Brunt-Vaisala  
    1288       !!      frequency at the time-step of the input arguments 
    1289       !! 
    1290       !! ** Method  :   call eos_ab to obtain alpha and beta coefficients 
    1291       !!                then interpolate them vertically, and compute pn2 
    1292       !!        Macro-tasked on horizontal slab (jk-loop) 
    1293       !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr 
    1294       !!      and is never used at this level. 
    1295       !! 
    1296       !! ** Action  : - pn2 : the brunt-vaisala frequency 
    1297       !! 
    1298       !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
    1299       !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    1300       !!                McDougall, J. Phys. Oceano., 1987 
    1301       !!---------------------------------------------------------------------- 
    1302       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    1303       !                                                               ! 2 : salinity               [psu] 
    1304       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2   ! Brunt-Vaisala frequency    [s-1] 
    1305       !! 
    1306       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    1307       REAL(wp) ::   zgde3w, za, zb, zr1          ! temporary scalars 
    1308 #if defined key_zdfddm 
    1309       REAL(wp) ::   zds   ! local scalars 
    1310 #endif 
    1311       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zab 
    1312       !!---------------------------------------------------------------------- 
    1313       ! 
    1314       ! 
    1315       IF( nn_timing == 1 ) CALL timing_start('eos_bn2_ab') 
    1316       ! 
    1317       ! pn2 : interior points only (2=< jk =< jpkm1 ) 
    1318       ! -------------------------- 
    1319       ! 
    1320       CALL wrk_alloc( jpi, jpj, jpk, jpts, zab ) 
    1321       ! 
    1322       CALL eos_ab( pts, zab ) 
    1323       ! 
    1324       DO jk = 2, jpkm1 
    1325          DO jj = 1, jpj 
    1326             DO ji = 1, jpi 
    1327                ! 
    1328                zgde3w = 1._wp / fse3w(ji,jj,jk) 
    1329                zr1 = ( fsdepw(ji,jj,jk) - fsdept(ji,jj,jk) ) * zgde3w 
    1330                ! 
    1331                za = ( zab(ji,jj,jk,jp_tem) +   & 
    1332                   &   ( zab(ji,jj,jk-1,jp_tem) - zab(ji,jj,jk,jp_tem) ) * zr1 )   & 
    1333                   &   * tmask(ji,jj,jk) 
    1334                zb = ( zab(ji,jj,jk,jp_sal) +   & 
    1335                   &   ( zab(ji,jj,jk-1,jp_sal) - zab(ji,jj,jk,jp_sal) ) * zr1 )   & 
    1336                   &   * tmask(ji,jj,jk) 
    1337                ! N2 
    1338                pn2(ji,jj,jk) = grav * zgde3w                    &   ! N^2 
    1339                   &          * (   za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    1340                   &              - zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
    1341 #if defined key_zdfddm 
    1342                zds = zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )      ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    1343                IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    1344                rrau(ji,jj,jk) = za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    1345 #endif 
    1346             END DO 
    1347          END DO 
    1348       END DO 
    1349       ! 
    1350  
    1351       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk ) 
    1352 #if defined key_zdfddm 
    1353       IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk ) 
    1354 #endif 
    1355       ! 
    1356       CALL wrk_dealloc( jpi, jpj, jpk, jpts, zab ) 
    1357       ! 
    1358       IF( nn_timing == 1 ) CALL timing_stop('eos_bn2_ab') 
    1359       ! 
    1360    END SUBROUTINE eos_bn2_ab 
    1361  
    1362  
    1363    SUBROUTINE eos_bn2( pts, pab, pn2 ) 
    1364       !!---------------------------------------------------------------------- 
    1365       !!                  ***  ROUTINE eos_bn2  *** 
    1366       !! 
    1367       !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the  
    1368       !!                time-step of the input arguments 
    1369       !! 
    1370       !! ** Method  :   interpolate alpha/beta vertically, and compute pn2 
    1371       !!        Macro-tasked on horizontal slab (jk-loop) 
    1372       !!      N.B. N^2 is set to zero at the first level (JK=1) in inidtr 
    1373       !!      and is never used at this level. 
    1374       !! 
    1375       !! ** Action  : - pn2 : the brunt-vaisala frequency 
    1376       !! 
    1377       !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
    1378       !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    1379       !!                McDougall, J. Phys. Oceano., 1987 
    1380       !!---------------------------------------------------------------------- 
    1381       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu] 
    1382       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! partial derivative of in situ  
    1383       !                                                              ! density with respect to T and S [1/Celcius,1/psu] 
    1384       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency    [1/s^2] 
    1385       !! 
    1386       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    1387       REAL(wp) ::   zgde3w, za, zb, zr1          ! temporary scalars 
    1388 #if defined key_zdfddm 
    1389       REAL(wp) ::   zds   ! local scalars 
    1390 #endif 
    1391       !!---------------------------------------------------------------------- 
    1392       ! 
    1393       ! 
    1394       IF( nn_timing == 1 ) CALL timing_start('eos_bn2') 
    1395       ! 
    1396       ! pn2 : interior points only (2=< jk =< jpkm1 ) 
    1397       ! -------------------------- 
    1398       ! 
    1399       DO jk = 2, jpkm1 
    1400          DO jj = 1, jpj 
    1401             DO ji = 1, jpi 
    1402                zgde3w = 1._wp / fse3w(ji,jj,jk) 
    1403                zr1 = ( fsdepw(ji,jj,jk) - fsdept(ji,jj,jk) ) * zgde3w 
    1404                ! 
    1405                za = ( pab(ji,jj,jk,jp_tem) + ( pab(ji,jj,jk-1,jp_tem) - pab(ji,jj,jk,jp_tem) ) * zr1 ) * tmask(ji,jj,jk) 
    1406                zb = ( pab(ji,jj,jk,jp_sal) + ( pab(ji,jj,jk-1,jp_sal) - pab(ji,jj,jk,jp_sal) ) * zr1 ) * tmask(ji,jj,jk) 
    1407                ! 
    1408                pn2(ji,jj,jk) = grav * zgde3w                    &   ! N^2 
    1409                   &          * (   za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    1410                   &              - zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
    1411 #if defined key_zdfddm 
    1412                zds = zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )      ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    1413                IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    1414                rrau(ji,jj,jk) = za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    1415 #endif 
    1416             END DO 
    1417          END DO 
    1418       END DO 
    1419       ! 
    1420  
    1421       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk ) 
    1422 #if defined key_zdfddm 
    1423       IF(ln_ctl)   CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk ) 
    1424 #endif 
    1425       ! 
    1426       IF( nn_timing == 1 )   CALL timing_stop('eos_bn2') 
    1427       ! 
    1428    END SUBROUTINE eos_bn2 
    1429  
    1430  
    1431  
    1432    FUNCTION tfreez( psal ) RESULT( ptf ) 
    1433       !!---------------------------------------------------------------------- 
    1434       !!                 ***  ROUTINE tfreez  *** 
    1435       !! 
    1436       !! ** Purpose :   Compute the sea surface freezing temperature [Celcius] 
    1437       !! 
    1438       !! ** Method  :   UNESCO freezing point at the surface (pressure = 0???) 
    1439       !!       freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p 
    1440       !!       checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars 
    1441       !! 
    1442       !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978 
    1443       !!---------------------------------------------------------------------- 
    1444       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu] 
    1445       ! Leave result array automatic rather than making explicitly allocated 
    1446       REAL(wp), DIMENSION(jpi,jpj)                ::   ptf    ! freezing temperature [Celcius] 
    1447       !!---------------------------------------------------------------------- 
    1448       ! 
    1449       ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
    1450          &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:) 
    1451       ! 
    1452    END FUNCTION tfreez 
    1453788 
    1454789 
     
    1495830         IF(lwp) WRITE(numout,*) '          use of simplified eos:    rhd(T,S,z) = ' 
    1496831         IF(lwp) WRITE(numout,*) '             -alph0*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta0*(S-S0) + gamm0*z' 
    1497          IF( lk_zdfddm .AND. rn_beta0==0. )   & 
    1498               &   CALL ctl_stop( '          double diffusive mixing parameterization requires that rn_beta<>0') 
    1499          offset = 1027._wp / rau0 - 1._wp 
     832         ! 
     833         offset = 1027._wp * r1_rau0 - 1._wp             ! offset value so that rho(10,35,0)=1027kg/m3 
     834         ! 
     835         rab_b(:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:)   ! beta = constant : set one for all 
     836         rab_n(:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:) 
     837         ! 
     838         IF(lwp) WRITE(numout,*) '          offset value so that rho(10,35,0)=1027kg/m3 :   offset = ', offset 
     839         ! 
     840         IF( rn_beta0 == 0._wp ) THEN      ! rho is not a function of salinty 
     841            IF( lk_ldfslp )   CALL ctl_stop( 'use of isoneutral mixing requires that rn_beta/=0') 
     842            IF( lk_zdfddm )   CALL ctl_stop( 'use of double diffusive mixing requires that rn_beta/=0') 
     843         ENDIF 
    1500844         ! 
    1501845      CASE DEFAULT                     !==  ERROR in nn_eos  ==! 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90

    r3876 r3894  
    127127      REAL(wp) ::   zupsut, zcenut, zupst                 !   -      - 
    128128      REAL(wp) ::   zupsvt, zcenvt, zcent, zice           !   -      - 
    129       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztfreez  
    130       REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind 
     129      REAL(wp), POINTER, DIMENSION(:,:)   :: zfzp         ! 2D workspace 
     130      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind    ! 3D workspace 
    131131      !!---------------------------------------------------------------------- 
    132132      ! 
    133133      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen2') 
    134134      ! 
    135       CALL wrk_alloc( jpi, jpj, ztfreez ) 
     135      CALL wrk_alloc( jpi, jpj, zfzp ) 
    136136      CALL wrk_alloc( jpi, jpj, jpk, zwz, zind ) 
    137137      ! 
     
    168168!!gm  not strickly exact : the freezing point should be computed at each ocean levels... 
    169169!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations 
    170       ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) ) 
     170      zfzp(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) 
    171171      DO jk = 1, jpk 
    172172         DO jj = 1, jpj 
    173173            DO ji = 1, jpi 
    174174               !                                        ! below ice covered area (if tn < "freezing"+0.1 ) 
    175                IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0 
    176                ELSE                                                      ;   zice = 0.e0 
     175               IF( tsn(ji,jj,jk,jp_tem) <= zfzp(ji,jj) + 0.1 ) THEN   ;   zice = 1._wp 
     176               ELSE                                                   ;   zice = 0._wp 
    177177               ENDIF 
    178178               zind(ji,jj,jk) = MAX (   & 
     
    280280      ENDIF 
    281281      ! 
    282       CALL wrk_dealloc( jpi, jpj, ztfreez ) 
     282      CALL wrk_dealloc( jpi, jpj, zfzp ) 
    283283      CALL wrk_dealloc( jpi, jpj, jpk, zwz, zind ) 
    284284      ! 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

    r3876 r3894  
    3030   USE lib_mpp        ! distribued memory computing 
    3131   USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    32 !!gm   USE trabbl         ! tracers: bottom boundary layer 
    33 !!gm   USE eosbn2          ! equation of state 
    3432 
    3533   IMPLICIT NONE 
    3634   PRIVATE 
    3735 
    38    PUBLIC   tra_adv_muscl  ! routine called by step.F90 
    39  
    40    LOGICAL  :: l_trd                        ! flag to compute trends 
    41    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits 
    42    !                                                             !  and in closed seas (orca 2 and 4 configurations) 
    43    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xind         !: mixed upstream/centered index 
     36   PUBLIC   tra_adv_muscl   ! routine called by step.F90 
     37 
     38   LOGICAL  ::   l_trd   ! flag to compute trends 
     39    
     40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits 
     41   !                                                           !  and in closed seas (orca 2 and 4 configurations) 
     42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index 
     43    
    4444   !! * Substitutions 
    4545#  include "domzgr_substitute.h90" 
     
    5252CONTAINS 
    5353 
    54    SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 
    55       &                                        ptb, pta, kjpt, ld_msc_ups ) 
     54   SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn,   & 
     55      &                                                ptb, pta, kjpt, ld_msc_ups ) 
    5656      !!---------------------------------------------------------------------- 
    5757      !!                    ***  ROUTINE tra_adv_muscl  *** 
     
    6969      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 
    7070      !!---------------------------------------------------------------------- 
    71       USE oce     , ONLY:   zwx   => ua    , zwy   => va          ! (ua,va) used as workspace 
     71      USE oce, ONLY:   zwx   => ua    , zwy   => va          ! (ua,va) used as workspace 
    7272      ! 
    7373      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
     
    8080      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb             ! before tracer field 
    8181      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
    82  
    83       ! 
    84       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     82      ! 
     83      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices 
     84      INTEGER  ::   ierr                      ! local integer 
    8585      REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars 
    8686      REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      - 
    8787      REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      - 
    88       REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy 
    89       INTEGER  ::   ierr 
     88      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zslpx, zslpy   ! workspace 
    9089      !!---------------------------------------------------------------------- 
    9190      ! 
     
    9493      CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy ) 
    9594      ! 
    96  
    9795      IF( kt == kit000 )  THEN 
    9896         IF(lwp) WRITE(numout,*) 
     
    121119 
    122120         ! 
    123          ! Upstream / centered scheme indicator 
     121         ! Upstream / MUSCL scheme indicator 
    124122         ! ------------------------------------ 
    125          xind(:,:,:) = 1._wp                             ! set equal to 1 where up-stream is not needed 
    126          ! 
    127123         IF( ld_msc_ups )  THEN 
    128124            DO jk = 1, jpk 
    129                DO jj = 1, jpj 
    130                   DO ji = 1, jpi 
    131                      xind(ji,jj,jk) = 1  - MAX (           & 
    132                         rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows) 
    133                         upsmsk(ji,jj) ) * tmask(ji,jj,jk)     ! some of some straits 
    134                   END DO 
    135                END DO 
     125               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed 
     126                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows) 
     127                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 near some straits 
    136128            END DO 
    137129         ENDIF  
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r3876 r3894  
    1212   !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC 
    1313   !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
     14   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta 
    1415   !!---------------------------------------------------------------------- 
    1516#if   defined key_trabbl   ||   defined key_esopa 
     
    2930   USE eosbn2         ! equation of state 
    3031   USE trd_oce        ! trends: ocean variables 
    31    USE trdtra         ! trends manager: tracers  
     32   USE trdtra         ! trends manager: tracers 
     33   ! 
    3234   USE iom            ! IOM server                
    3335   USE in_out_manager ! I/O manager 
     
    3638   USE wrk_nemo       ! Memory Allocation 
    3739   USE timing         ! Timing 
    38  
     40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3941 
    4042   IMPLICIT NONE 
     
    4951   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag 
    5052 
    51    !                                           !!* Namelist nambbl * 
    52    INTEGER , PUBLIC ::   nn_bbl_ldf = 0         !: =1   : diffusive bbl or not (=0) 
    53    INTEGER , PUBLIC ::   nn_bbl_adv = 0         !: =1/2 : advective bbl or not (=0) 
    54    !                                            !  =1 : advective bbl using the bottom ocean velocity 
    55    !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho) 
    56    REAL(wp), PUBLIC ::   rn_ahtbbl  = 1.e3_wp   !: along slope bbl diffusive coefficient [m2/s] 
    57    REAL(wp), PUBLIC ::   rn_gambbl  = 10.0_wp   !: lateral coeff. for bottom boundary layer scheme [s] 
    58  
    59    LOGICAL , PUBLIC ::   l_bbl                  !: flag to compute bbl diffu. flux coef and transport 
     53   !                                             !!* Namelist nambbl * 
     54   INTEGER , PUBLIC ::   nn_bbl_ldf = 0           !: =1   : diffusive bbl or not (=0) 
     55   INTEGER , PUBLIC ::   nn_bbl_adv = 0           !: =1/2 : advective bbl or not (=0) 
     56   !                                              !  =1 : advective bbl using the bottom ocean velocity 
     57   !                                              !  =2 :     -      -  using utr_bbl proportional to grad(rho) 
     58   REAL(wp), PUBLIC ::   rn_ahtbbl  = 1.e3_wp     !: along slope bbl diffusive coefficient [m2/s] 
     59   REAL(wp), PUBLIC ::   rn_gambbl  = 10.0_wp     !: lateral coeff. for bottom boundary layer scheme [s] 
     60 
     61   LOGICAL , PUBLIC ::   l_bbl                    !: flag to compute bbl diffu. flux coef. and transport 
    6062 
    6163   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer 
     
    105107      !!---------------------------------------------------------------------- 
    106108      INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
    107       !! 
     109      ! 
    108110      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    109111      !!---------------------------------------------------------------------- 
     
    111113      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl') 
    112114      ! 
    113       IF( l_trdtra )   THEN                        !* Save ta and sa trends 
     115      IF( l_trdtra )   THEN                         !* Save ta and sa trends 
    114116         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    115117         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
     
    117119      ENDIF 
    118120 
    119       IF( l_bbl )  CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
    120  
    121       IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl 
     121      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
     122 
     123      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl 
    122124         ! 
    123125         CALL tra_bbl_dif( tsb, tsa, jpts ) 
    124126         IF( ln_ctl )  & 
    125127         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, & 
    126          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     128            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    127129         ! lateral boundary conditions ; just need for outputs 
    128130         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. ) 
     
    132134      END IF 
    133135 
    134       IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl 
     136      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl 
    135137         ! 
    136138         CALL tra_bbl_adv( tsb, tsa, jpts ) 
    137139         IF(ln_ctl)   & 
    138140         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   & 
    139          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     141            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    140142         ! lateral boundary conditions ; just need for outputs 
    141143         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. ) 
     
    165167      !!                advection terms. 
    166168      !! 
    167       !! ** Method  : 
    168       !!        * diffusive bbl (nn_bbl_ldf=1) : 
     169      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) : 
    169170      !!        When the product grad( rho) * grad(h) < 0 (where grad is an 
    170171      !!      along bottom slope gradient) an additional lateral 2nd order 
     
    196197      DO jn = 1, kjpt                                     ! tracer loop 
    197198         !                                                ! =========== 
    198 #  if defined key_vectopt_loop 
    199          DO jj = 1, 1   ! vector opt. (forced unrolling) 
    200             DO ji = 1, jpij 
    201 #else 
    202199         DO jj = 1, jpj 
    203200            DO ji = 1, jpi 
    204 #endif 
    205                ik = mbkt(ji,jj)                        ! bottom T-level index 
    206                zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S 
     201               zptb(ji,jj) = ptb(ji,jj,mbkt(ji,jj),jn)       ! bottom before T and S 
    207202            END DO 
    208203         END DO 
    209          !                                                ! Compute the trend 
    210 #  if defined key_vectopt_loop 
    211          DO jj = 1, 1   ! vector opt. (forced unrolling) 
    212             DO ji = jpi+1, jpij-jpi-1 
    213 #  else 
    214          DO jj = 2, jpjm1 
     204         !                
     205         DO jj = 2, jpjm1                                    ! Compute the trend 
    215206            DO ji = 2, jpim1 
    216 #  endif 
    217                ik = mbkt(ji,jj)                            ! bottom T-level index 
     207               ik = mbkt(ji,jj)                              ! bottom T-level index 
    218208               zbtr = r1_e1e2t(ji,jj)  / fse3t(ji,jj,ik) 
    219209               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         & 
     
    264254      DO jn = 1, kjpt                                            ! tracer loop 
    265255         !                                                       ! =========== 
    266 # if defined key_vectopt_loop 
    267          DO jj = 1, 1 
    268             DO ji = 1, jpij-jpi-1   ! vector opt. (forced unrolling) 
    269 # else 
    270256         DO jj = 1, jpjm1 
    271257            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west 
    272 # endif 
    273258               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection 
    274259                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf) 
     
    333318      !!                advection terms. 
    334319      !! 
    335       !! ** Method  : 
    336       !!        * diffusive bbl (nn_bbl_ldf=1) : 
     320      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) : 
    337321      !!        When the product grad( rho) * grad(h) < 0 (where grad is an 
    338322      !!      along bottom slope gradient) an additional lateral 2nd order 
     
    342326      !!      a downslope velocity of 20 cm/s if the condition for slope 
    343327      !!      convection is satified) 
    344       !!        * advective bbl (nn_bbl_adv=1 or 2) : 
     328      !!              * advective bbl (nn_bbl_adv=1 or 2) : 
    345329      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity 
    346330      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation 
     
    353337      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 
    354338      !!---------------------------------------------------------------------- 
    355       ! 
    356339      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index 
    357       INTEGER         , INTENT(in   ) ::   kit000          ! first time step index 
     340      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index 
    358341      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    359       !! 
    360       INTEGER  ::   ji, jj                    ! dummy loop indices 
    361       INTEGER  ::   ik                        ! local integers 
    362       INTEGER  ::   iis , iid , ijs , ijd     !   -       - 
    363       INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       - 
    364       REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars 
    365       REAL(wp) ::   zgdrho, zt, zs, zh        !   -      - 
    366       !! 
    367       REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function 
    368       REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep 
    369       !!----------------------- zv_bbl----------------------------------------------- 
    370       ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients 
    371       ! ================            pft :  potential temperature in degrees celcius 
    372       !                             pfs :  salinity anomaly (s-35) in psu 
    373       !                             pfh :  depth in meters 
    374       ! nn_eos = 0  (Jackett and McDougall 1994 formulation) 
    375       fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta 
    376          ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    & 
    377                                    - 0.203814e-03 ) * pft                    & 
    378                                    + 0.170907e-01 ) * pft                    & 
    379                                    + 0.665157e-01                            & 
    380          +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   & 
    381          +  ( ( - 0.302285e-13 * pfh                                         & 
    382                 - 0.251520e-11 * pfs                                         & 
    383                 + 0.512857e-12 * pft * pft          ) * pfh                  & 
    384                                      - 0.164759e-06   * pfs                  & 
    385              +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  & 
    386                                      + 0.380374e-04 ) * pfh 
    387       fsbeta( pft, pfs, pfh ) =                                              &   ! beta 
    388          ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      & 
    389                                  - 0.301985e-05 ) * pft                      & 
    390                                  + 0.785567e-03                              & 
    391          + (     0.515032e-08 * pfs                                          & 
    392                + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   & 
    393                +(  (   0.121551e-17 * pfh                                    & 
    394                      - 0.602281e-15 * pfs                                    & 
    395                      - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             & 
    396                                           + 0.408195e-10   * pfs             & 
    397                  + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             & 
    398                                           - 0.121555e-07 ) * pfh 
     342      ! 
     343      INTEGER  ::   ji, jj                      ! dummy loop indices 
     344      INTEGER  ::   ik, ikip1, ikjp1            ! local integers 
     345      INTEGER  ::   iis, iid, ikus, ikud        !   -       - 
     346      INTEGER  ::   ijs, ijd, ikvs, ikvd        !   -       - 
     347      REAL(wp) ::   za, zt, zsign , zgbbl       ! local scalars 
     348      REAL(wp) ::   zb, zs, zsigna, zgdrho, zh  !   -      - 
     349      REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts, zab   ! 3D workspace 
     350      REAL(wp), DIMENSION(jpi,jpj)      :: zub, zvb   ! 2D workspace 
    399351      !!---------------------------------------------------------------------- 
    400352      ! 
    401353      IF( nn_timing == 1 )  CALL timing_start( 'bbl') 
    402       ! 
    403       CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
    404354      ! 
    405355      IF( kt == kit000 )  THEN 
     
    408358         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~' 
    409359      ENDIF 
    410       !                                        !* bottom temperature, salinity, velocity and depth 
    411 #if defined key_vectopt_loop 
    412       DO jj = 1, 1   ! vector opt. (forced unrolling) 
    413          DO ji = 1, jpij 
    414 #else 
     360      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity) 
    415361      DO jj = 1, jpj 
    416362         DO ji = 1, jpi 
    417 #endif 
    418             ik = mbkt(ji,jj)                        ! bottom T-level index 
    419             ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S 
    420             zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1) 
    421             zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth 
     363            ik = mbkt(ji,jj)                             ! bottom T-level index 
     364            zts(ji,jj,jp_tem) = tsb  (ji,jj,ik,jp_tem)   ! bottom before T and S 
     365            zts(ji,jj,jp_sal) = tsb  (ji,jj,ik,jp_sal) 
     366            zab(ji,jj,jp_tem) = rab_b(ji,jj,ik,jp_tem)   ! bottom alpha and beta 
     367            zab(ji,jj,jp_sal) = rab_b(ji,jj,ik,jp_sal) 
    422368            ! 
    423             zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity 
     369            zub(ji,jj) = un(ji,jj,mbku(ji,jj))           ! bottom velocity 
    424370            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
    425371         END DO 
    426372      END DO 
    427  
     373      ! 
    428374      !                                   !-------------------! 
    429375      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   ! 
    430376         !                                !-------------------! 
    431377         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 
    432             DO ji = 1, jpim1 
    433                !                                                ! i-direction 
    434                zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth 
    435                zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    436                zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    437                !                                                         ! masked bbl i-gradient of density 
    438                zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
    439                   &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1) 
     378            DO ji = 1, fs_jpim1   ! vector opt. 
     379               !                                                   ! i-direction 
     380               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point 
     381               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     382               !                                                         ! 2*masked bottom density gradient 
     383               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
     384                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
    440385               ! 
    441                zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope ) 
    442                ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff. 
     386               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope ) 
     387               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff. 
    443388               ! 
    444                !                                                ! j-direction 
    445                zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth 
    446                zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 
    447                zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 
    448                !                                                         ! masked bbl j-gradient of density 
    449                zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
    450                   &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1) 
     389               !                                                   ! j-direction 
     390               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point 
     391               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     392               !                                                         ! 2*masked bottom density gradient 
     393               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    & 
     394                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1) 
    451395               ! 
    452                zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope ) 
     396               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope ) 
    453397               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj) 
    454                ! 
    455398            END DO 
    456399         END DO 
     
    466409            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0 
    467410               DO ji = 1, fs_jpim1   ! vector opt. 
    468                   !                                               ! i-direction 
    469                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth 
    470                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    471                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    472                   !                                                           ! masked bbl i-gradient of density 
    473                   zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
    474                      &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1) 
    475                   ! 
    476                   zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope 
    477                   zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope 
    478                   ! 
    479                   !                                                           ! bbl velocity 
     411                  !                                                  ! i-direction 
     412                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
     413                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     414                  !                                                          ! 2*masked bottom density gradient  
     415                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
     416                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
     417                  ! 
     418                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope 
     419                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope 
     420                  ! 
     421                  !                                                          ! bbl velocity 
    480422                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj) 
    481423                  ! 
    482                   !                                               ! j-direction 
    483                   zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth 
    484                   zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 
    485                   zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 
    486                   !                                                           ! masked bbl j-gradient of density 
    487                   zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
    488                      &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1) 
    489                   zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope 
    490                   zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope 
    491                   ! 
    492                   !                                                           ! bbl velocity 
     424                  !                                                  ! j-direction 
     425                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point 
     426                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     427                  !                                                          ! 2*masked bottom density gradient 
     428                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    & 
     429                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1) 
     430                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope 
     431                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope 
     432                  ! 
     433                  !                                                          ! bbl transport 
    493434                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj) 
    494435               END DO 
     
    499440            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down 
    500441               DO ji = 1, fs_jpim1   ! vector opt. 
    501                   !                                         ! i-direction 
     442                  !                                                  ! i-direction 
    502443                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf) 
    503444                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 
    504445                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj) 
    505446                  ! 
    506                   !                                             ! mid-depth density anomalie (up-slope minus down-slope) 
    507                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth 
    508                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    509                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    510                   zgdrho =    fsbeta( zt, zs, zh )                                    & 
    511                      &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    & 
    512                      &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1) 
    513                   zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep 
    514                   ! 
    515                   !                                             ! bbl transport (down-slope direction) 
     447                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
     448                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     449                  !                                                          !   masked bottom density gradient 
     450                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    & 
     451                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1) 
     452                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep 
     453                  ! 
     454                  !                                                          ! bbl transport (down-slope direction) 
    516455                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) ) 
    517456                  ! 
    518                   !                                         ! j-direction 
     457                  !                                                  ! j-direction 
    519458                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf) 
    520459                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 
    521460                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj) 
    522461                  ! 
    523                   !                                             ! mid-depth density anomalie (up-slope minus down-slope) 
    524                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth 
    525                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0 
    526                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) ) 
    527                   zgdrho =    fsbeta( zt, zs, zh )                                    & 
    528                      &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    & 
    529                      &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1) 
    530                   zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep 
    531                   ! 
    532                   !                                             ! bbl transport (down-slope direction) 
     462                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point 
     463                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     464                  !                                                          !   masked bottom density gradient 
     465                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    & 
     466                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1) 
     467                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep 
     468                  ! 
     469                  !                                                          ! bbl transport (down-slope direction) 
    533470                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) ) 
    534471               END DO 
     
    537474         ! 
    538475      ENDIF 
    539       ! 
    540       CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
    541476      ! 
    542477      IF( nn_timing == 1 )  CALL timing_stop( 'bbl') 
     
    606541 
    607542                                        !* sign of grad(H) at u- and v-points 
    608       mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0. 
     543      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0 
    609544      DO jj = 1, jpjm1 
    610545         DO ji = 1, jpim1 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

    r3876 r3894  
    8181         inpci = 0 
    8282 
    83          CALL eos( tsa, rhd, zrhop )         ! Potential density 
     83         CALL eos( tsa, fsdept(:,:,:), rhd, zrhop )         ! Potential density 
    8484 
    8585         IF( l_trdtra )   THEN                    !* Save ta and sa trends 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r3294 r3894  
    7474      !!          Idem for di(s) and dj(s)           
    7575      !! 
    76       !!      For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at  
    77       !!      the good depth zh from interpolated T and S for the different 
    78       !!      formulation of the equation of state (eos). 
     76      !!      For rho, we call eos which will compute rd~(t~,s~) at the right 
     77      !!      depth zh from interpolated T and S for the different formulations 
     78      !!      of the equation of state (eos). 
    7979      !!      Gradient formulation for rho : 
    80       !!          di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 
     80      !!          di(rho) = rd~ - rd(i,j,k)   or  rd(i+1,j,k) - rd~ 
    8181      !! 
    8282      !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
    8383      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points  
    8484      !!---------------------------------------------------------------------- 
    85       ! 
    8685      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    8786      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers 
    8887      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
    8988      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    90       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    91       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad. of prd at u- & v-pts  
     89      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
     90      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad. of prd at u- & v-pts  
    9291      ! 
    9392      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    9493      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
    9594      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
    96       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zri, zrj, zhi, zhj 
    97       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zti, ztj    ! interpolated value of tracer 
     95      REAL(wp), DIMENSION(jpi,jpj,1)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     96      REAL(wp), DIMENSION(jpi,jpj,1,kjpt) ::  zti, ztj             !  
    9897      !!---------------------------------------------------------------------- 
    9998      ! 
    10099      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
    101       ! 
    102       CALL wrk_alloc( jpi, jpj,       zri, zrj, zhi, zhj )  
    103       CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj           )  
    104100      ! 
    105101      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
     
    121117                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    122118                  ! interpolated values of tracers 
    123                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     119                  zti(ji,jj,1,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
    124120                  ! gradient of  tracers 
    125                   pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     121                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,1,jn) - pta(ji,jj,iku,jn) ) 
    126122               ELSE                           ! case 2 
    127123                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    128124                  ! interpolated values of tracers 
    129                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     125                  zti(ji,jj,1,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    130126                  ! gradient of tracers 
    131                   pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     127                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,1,jn) ) 
    132128               ENDIF 
    133129               ! 
     
    136132                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    137133                  ! interpolated values of tracers 
    138                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     134                  ztj(ji,jj,1,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    139135                  ! gradient of tracers 
    140                   pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     136                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,1,jn) - pta(ji,jj,ikv,jn) ) 
    141137               ELSE                           ! case 2 
    142138                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    143139                  ! interpolated values of tracers 
    144                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     140                  ztj(ji,jj,1,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    145141                  ! gradient of tracers 
    146                   pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     142                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,1,jn) ) 
    147143               ENDIF 
    148144# if ! defined key_vectopt_loop 
     
    167163               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    168164               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    169                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
    170                ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
    171                ENDIF 
    172                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
    173                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
     165               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj,1) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
     166               ELSE                        ;   zhi(ji,jj,1) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
     167               ENDIF 
     168               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj,1) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
     169               ELSE                        ;   zhj(ji,jj,1) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
    174170               ENDIF 
    175171# if ! defined key_vectopt_loop 
     
    195191               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    196192               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    197                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    198                ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
    199                ENDIF 
    200                IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
    201                ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
     193               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj, 1 ) - prd(ji,jj,iku) )   ! i: 1 
     194               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj, 1 ) )   ! i: 2 
     195               ENDIF 
     196               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj  , 1 ) - prd(ji,jj,ikv) )   ! j: 1 
     197               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj, 1 ) )   ! j: 2 
    202198               ENDIF 
    203199# if ! defined key_vectopt_loop 
     
    209205      END IF 
    210206      ! 
    211       CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj )  
    212       CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj           )  
    213       ! 
    214207      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
    215208      ! 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90

    r3876 r3894  
    212212         zkepe(:,:,:) = 0._wp 
    213213    
    214          CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities 
     214         CALL eos( tsn, fsdept(:,:,:), rhd, rhop )       ! now potential and in situ densities 
    215215 
    216216         zcof = 0.5_wp / rau0             ! Density flux at w-point 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r3610 r3894  
    66   !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing 
    77   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
    8    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     8   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     9   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_zdfddm   ||   defined key_esopa 
     
    1819   USE dom_oce         ! ocean space and time domain variables  
    1920   USE zdf_oce         ! ocean vertical physics variables 
     21   ! 
    2022   USE in_out_manager  ! I/O manager 
    2123   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    3436   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag 
    3537 
    36    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs    !: salinity vertical diffusivity coeff. at w-point 
    37    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   rrau   !: heat/salt buoyancy flux ratio 
     38   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs   !: salinity vertical diffusivity coeff. at w-point 
    3839 
    3940   !                                  !!* Namelist namzdf_ddm : double diffusive mixing * 
     
    4243 
    4344   !! * Substitutions 
     45#  include "domzgr_substitute.h90" 
    4446#  include "vectopt_loop_substitute.h90" 
    4547   !!---------------------------------------------------------------------- 
    46    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     48   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
    4749   !! $Id$ 
    4850   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    5456      !!                ***  ROUTINE zdf_ddm_alloc  *** 
    5557      !!---------------------------------------------------------------------- 
    56       ALLOCATE( avs(jpi,jpj,jpk), rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc ) 
     58      ALLOCATE( avs(jpi,jpj,jpk) , STAT= zdf_ddm_alloc ) 
    5759      ! 
    5860      IF( lk_mpp             )   CALL mpp_sum ( zdf_ddm_alloc ) 
     
    7173      !!      diffusive mixing (i.e. salt fingering and diffusive layering) 
    7274      !!      following Merryfield et al. (1999). The rate of double diffusive  
    73       !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S] 
    74       !!      which is computed in rn2.F 
     75      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]): 
    7576      !!         * salt fingering (Schmitt 1981): 
    76       !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 ) 
    77       !!      for Rrau > 1 and rn2 > 0 : zavfs = O 
    78       !!      otherwise                : zavft = 0.7 zavs / Rrau 
     77      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 ) 
     78      !!      for R > 1 and rn2 > 0 : zavfs = O 
     79      !!      otherwise                : zavft = 0.7 zavs / R 
    7980      !!         * diffusive layering (Federov 1988): 
    80       !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6   
    81       !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) ) 
     81      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) ) 
    8282      !!      otherwise                   : zavdt = 0  
    83       !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85) 
    84       !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau       
     83      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85) 
     84      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R       
    8585      !!      otherwise                     : zavds = 0  
    8686      !!         * update the eddy diffusivity: 
     
    9696      ! 
    9797      INTEGER  ::   ji, jj , jk     ! dummy loop indices 
    98       REAL(wp) ::   zinr, zrr       ! temporary scalars 
    99       REAL(wp) ::   zavft, zavfs    !    -         - 
    100       REAL(wp) ::   zavdt, zavds    !    -         - 
    101       REAL(wp), POINTER, DIMENSION(:,:) ::   zmsks, zmskf, zmskd1, zmskd2, zmskd3 
     98      REAL(wp) ::   zaw, zbw, zrw   ! local scalars 
     99      REAL(wp) ::   zdt, zds 
     100      REAL(wp) ::   zinr, zrr       !   -      - 
     101      REAL(wp) ::   zavft, zavfs    !   -      - 
     102      REAL(wp) ::   zavdt, zavds    !   -      - 
     103      REAL(wp), POINTER, DIMENSION(:,:) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 
    102104      !!---------------------------------------------------------------------- 
    103105      ! 
    104106      IF( nn_timing == 1 )  CALL timing_start('zdf_ddm') 
    105107      ! 
    106       CALL wrk_alloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 
     108      CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 
    107109 
    108110      !                                                ! =============== 
     
    111113         ! Define the mask  
    112114         ! --------------- 
    113          rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau 
     115         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 
     116            DO ji = 1, jpi 
     117               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   & 
     118                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )  
     119               ! 
     120               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  ) * tmask(ji,jj,jk) 
     121               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  ) * tmask(ji,jj,jk) 
     122               ! 
     123               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) 
     124               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) )  
     125               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp 
     126               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau 
     127            END DO 
     128         END DO 
    114129 
    115130         DO jj = 1, jpj                                     ! indicators: 
     
    119134               ELSE                                       ;   zmsks(ji,jj) = 1._wp 
    120135               ENDIF 
    121                ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere             
    122                IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp 
     136               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere             
     137               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp 
    123138               ELSE                                       ;   zmskf(ji,jj) = 1._wp 
    124139               ENDIF 
    125140               ! diffusive layering indicators:  
    126                !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere 
    127                IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp 
     141               !     ! mskdl1=1 if 0< R <1; 0 elsewhere 
     142               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp 
    128143               ELSE                                       ;   zmskd1(ji,jj) = 1._wp 
    129144               ENDIF 
    130                !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere 
    131                IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp 
     145               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere 
     146               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp 
    132147               ELSE                                       ;   zmskd2(ji,jj) = 1._wp 
    133148               ENDIF 
    134                !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere 
    135                IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp 
    136                ELSE                                                         ;   zmskd3(ji,jj) = 1._wp 
     149               !   mskdl3=1 if 0.5< R <1; 0 elsewhere 
     150               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp 
     151               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp 
    137152               ENDIF 
    138153            END DO 
     
    149164!CDIR NOVERRCHK 
    150165            DO ji = 1, jpi 
    151                zinr = 1./rrau(ji,jj,jk) 
     166               zinr = 1._wp / zrau(ji,jj) 
    152167               ! salt fingering 
    153                zrr = rrau(ji,jj,jk)/rn_hsbfr 
     168               zrr = zrau(ji,jj) / rn_hsbfr 
    154169               zrr = zrr * zrr 
    155170               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj) 
     
    157172               ! diffusive layering 
    158173               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj) 
    159                zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   & 
    160                   &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  ) 
     174               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   & 
     175                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  ) 
    161176               ! add to the eddy viscosity coef. previously computed 
    162177               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90

    r3792 r3894  
    2626   USE phycst         ! physical constants 
    2727   USE eosbn2         ! equation of state 
    28    USE zdfddm         ! double diffusion mixing 
     28   USE zdfddm         ! double diffusion mixing (avs array) 
    2929   USE in_out_manager ! I/O manager 
    3030   USE lib_mpp        ! MPP library 
     
    3232   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3333   USE prtctl         ! Print control 
    34    USE trdmod_oce     ! ocean trends definition 
     34   USE trd_oce        ! ocean trends definition 
    3535   USE trdtra         ! tracers trends 
    3636   USE timing         ! Timing 
     
    246246      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t 
    247247#if defined key_zdfddm 
    248       REAL(wp) ::   zrrau, zds, zavdds, zavddt,zinr   ! double diffusion mixing 
    249       REAL(wp), POINTER, DIMENSION(:,:)   ::     zdifs 
     248      REAL(wp) ::   zrw, zkm1s                    ! local scalars 
     249      REAL(wp) ::   zrrau, zdt, zds, zavdds, zavddt, zinr   ! double diffusion mixing 
     250      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdifs 
    250251      REAL(wp), POINTER, DIMENSION(:)     ::   za2s, za3s, zkmps 
    251       REAL(wp) ::                            zkm1s 
    252252      REAL(wp), POINTER, DIMENSION(:,:)   ::   zblcs 
    253253      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdiffus 
     
    332332                  avt (ji,jj,jk) =  avt (ji,jj,jk) + rn_difri * zfri     
    333333               ENDIF 
     334               ! 
    334335#if defined key_zdfddm  
    335                avs (ji,jj,jk) =  avt (ji,jj,jk)               
    336                !  Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90 
    337                ! ------------------------------------------------------------------ 
    338                ! only retains positive value of rrau 
    339                zrrau = MAX( rrau(ji,jj,jk), epsln ) 
    340                zds   = tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) 
    341                IF( zrrau > 1. .AND. zds > 0.) THEN 
    342                   ! 
    343                   ! Salt fingering case. 
    344                   !--------------------- 
    345                   ! Compute interior diffusivity for double diffusive mixing of 
    346                   ! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001). 
    347                   ! After that set interior diffusivity for double diffusive mixing 
    348                   ! of temperature 
     336               ! 
     337               !  Double diffusion mixing ; NOT in ROUTINE zdfddm.F90 
     338               ! ------------------------- 
     339               avs (ji,jj,jk) =  avt (ji,jj,jk)    
     340 
     341               ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 
     342               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   & 
     343                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )  
     344               ! 
     345               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  ) * tmask(ji,jj,jk) 
     346               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  ) * tmask(ji,jj,jk) 
     347               ! 
     348               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) 
     349               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) )  
     350               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp 
     351               zrrau = MAX(  epsln , zdt / zds  )    ! only retains positive value of zrau 
     352               ! 
     353               IF( zrrau > 1. .AND. zds > 0.) THEN                        ! Salt fingering case. 
     354                  !                                                       !--------------------- 
     355                  ! Compute interior diffusivity for double diffusive mixing of salinity.  
     356                  ! Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001). 
     357                  ! After that set interior diffusivity for double diffusive mixing of temperature 
    349358                  zavdds = MIN( zrrau, Rrho0 ) 
    350359                  zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 ) 
     
    353362                  zavdds = difssf * zavdds  
    354363                  zavddt = 0.7 * zavdds 
    355                ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN 
    356364                  ! 
    357                   ! Diffusive convection case. 
    358                   !--------------------------- 
    359                   ! Compute interior diffusivity for double diffusive mixing of 
    360                   ! temperature (Marmorino and Caldwell, 1976);  
     365               ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN   ! Diffusive convection case. 
     366                  !                                                       !--------------------------- 
     367                  ! Compute interior diffusivity for double diffusive mixing of temperature (Marmorino and Caldwell, 1976);  
    361368                  ! Compute interior diffusivity for double diffusive mixing of salinity  
    362369                  zinr   = 1. / zrrau 
    363370                  zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) )  
    364371                  zavddt = difsdc * zavddt 
    365                   IF( zrrau < 0.5) THEN 
    366                      zavdds = zavddt * 0.15 * zrrau 
    367                   ELSE 
    368                      zavdds = zavddt * (1.85 * zrrau - 0.85 )  
     372                  IF( zrrau < 0.5) THEN   ;   zavdds = zavddt * 0.15 * zrrau 
     373                  ELSE                    ;   zavdds = zavddt * (1.85 * zrrau - 0.85 )  
    369374                  ENDIF 
    370375               ELSE 
     
    385390      !--------------------------------------------------------------------- 
    386391      DO jj = 2, jpjm1 
    387          DO ji = fs_2, fs_jpim1      
    388             IF( nn_eos < 1) THEN    
    389                zt     = tsn(ji,jj,1,jp_tem) 
    390                zs     = tsn(ji,jj,1,jp_sal) - 35.0 
    391                zh     = fsdept(ji,jj,1) 
    392                !  potential volumic mass 
    393                zrhos  = rhop(ji,jj,1) 
    394                zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta 
    395                   &                               - 0.203814e-03 ) * zt   & 
    396                   &                               + 0.170907e-01 ) * zt   & 
    397                   &   + 0.665157e-01                                      & 
    398                   &   +     ( - 0.678662e-05 * zs                         & 
    399                   &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   & 
    400                   &   +   ( ( - 0.302285e-13 * zh                         & 
    401                   &           - 0.251520e-11 * zs                         & 
    402                   &           + 0.512857e-12 * zt * zt           ) * zh   & 
    403                   &           - 0.164759e-06 * zs                         & 
    404                   &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   & 
    405                   &                               + 0.380374e-04 ) * zh 
    406  
    407                zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta 
    408                   &                            - 0.301985e-05 ) * zt      & 
    409                   &   + 0.785567e-03                                      & 
    410                   &   + (     0.515032e-08 * zs                           & 
    411                   &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     & 
    412                   &   +(  (   0.121551e-17 * zh                           & 
    413                   &         - 0.602281e-15 * zs                           & 
    414                   &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     & 
    415                   &                             + 0.408195e-10   * zs     & 
    416                   &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     & 
    417                   &                             - 0.121555e-07 ) * zh 
    418  
    419                zthermal = zbeta * zalbet / ( rcp * zrhos + epsln ) 
    420                zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs 
    421             ELSE 
    422                zrhos    = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) ) 
    423                zthermal = rn_alpha / ( rcp * zrhos + epsln ) 
    424                zhalin   = rn_beta * tsn(ji,jj,1,jp_sal) * rcs 
    425                zbeta    = rn_beta 
    426             ENDIF 
     392         DO ji = fs_2, fs_jpim1            
     393            zrhos    = rau0 * ( 1._wp + rhd(ji,jj,1) ) * tmask(ji,jj,1) 
     394            zthermal = rab_n(ji,jj,1,jp_tem) / ( rcp * zrhos + epsln ) 
     395            zbeta    = rab_n(ji,jj,1,jp_sal) 
     396            zhalin   = zbeta * tsn(ji,jj,1,jp_sal) * rcs 
     397            ! 
    427398            ! Radiative surface buoyancy force 
    428399            zBosol(ji,jj) = grav * zthermal * qsr(ji,jj) 
     
    435406            ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)                          & 
    436407               &             + sfx(ji,jj)                                     ) * rcs * tmask(ji,jj,1)  
    437          ENDDO 
    438       ENDDO 
     408         END DO 
     409      END DO 
    439410 
    440411      zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. )  
     
    447418            ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
    448419            zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos +  epsln ) ) 
    449          ENDDO 
    450       ENDDO 
     420         END DO 
     421      END DO 
    451422 
    452423!CDIR NOVERRCHK   
     
    12701241         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    12711242!!bug gm jpttdzdf ==> jpttkpp 
    1272          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdf, ztrdt ) 
    1273          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdf, ztrds ) 
     1243         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
     1244         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    12741245         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds ) 
    12751246      ENDIF 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r3625 r3894  
    435435 
    436436         ztpc = 0.e0 
    437          zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 
     437         zpc(:,:,:) = MAX( rn_n2min , rn2(:,:,:) ) * zav_tide(:,:,:) 
    438438         DO jk= 2, jpkm1 
    439439            DO jj = 1, jpj 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r3893 r3894  
    2929   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2     [s-2] 
    3030   ! 
    31    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units] 
     31   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: density anomaly rhd=(rho-rau0)/rau0  [no units] 
    3232   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhop   !: potential volumic mass                           [kg/m3] 
    3333 
  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/step.F90

    r3893 r3894  
    3838   !! * Substitutions 
    3939#  include "domzgr_substitute.h90" 
    40 #  include "zdfddm_substitute.h90" 
    41    !!---------------------------------------------------------------------- 
    42    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     40   !!---------------------------------------------------------------------- 
     41   !! NEMO/OPA 3.5 , NEMO Consortium (2010) 
    4342   !! $Id$ 
    4443   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    104103      ! Ocean physics update                (ua, va, tsa used as workspace) 
    105104      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    106                          CALL eos( tsb, rab_b )              ! before thermal/haline expansion coef.  
    107                          CALL bn2( tsb, rab_b, rn2b )        ! before Brunt-Vaisala frequency 
    108                          CALL eos( tsn, rab_n )              ! now thermal/haline expansion coef.  
    109                          CALL bn2( tsn, rab_n, rn2  )        ! now    Brunt-Vaisala frequency 
    110                           
     105      !                                               ! thermal/haline expansion coef. & Brunt-Vaisala frequency 
     106                         CALL eos( tsb, rab_b, rn2b )    ! before  
     107                         CALL eos( tsn, rab_n, rn2  )    ! now 
    111108      ! 
    112109      !  VERTICAL PHYSICS 
     
    142139      ! 
    143140      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing 
    144                          CALL eos( tsb, rhd )                ! before in situ density 
     141                         CALL eos( tsb, fsdept(:,:,:), rhd )              ! before in situ density 
    145142         IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
    146143            &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
     
    205202         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    206203                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
    207                              CALL eos    ( tsa, rhd, rhop )      ! Time-filtered in situ density for hpg computation 
     204                             CALL eos    ( tsa, fsdept(:,:,:), rhd, rhop )    ! Time-filtered in situ density for hpg computation 
    208205         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! zps: time filtered hor. derivative 
    209206            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
    210207 
    211208      ELSE                                                  ! centered hpg  (eos then time stepping) 
    212                              CALL eos    ( tsn, rhd, rhop )      ! now in situ density for hpg computation 
     209                             CALL eos    ( tsn, fsdept(:,:,:), rhd, rhop )    ! now in situ density for hpg computation 
    213210         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative 
    214211            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
Note: See TracChangeset for help on using the changeset viewer.