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

Changeset 3876


Ignore:
Timestamp:
2013-04-19T08:48:21+02:00 (11 years ago)
Author:
gm
Message:

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

Location:
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM
Files:
9 added
7 deleted
42 edited

Legend:

Unmodified
Added
Removed

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

    r3795 r3876  
    848848!!====================================================================== 
    849849!!   namnc4       netcdf4 chunking and compression settings             ("key_netcdf4") 
    850 !!   namtrd       dynamics and/or tracer trends                         ("key_trddyn","key_trdtra","key_trdmld") 
     850!!   namtrd       dynamics and/or tracer trends                          
    851851!!   namflo       float parameters                                      ("key_float") 
    852852!!   namptr       Poleward Transport Diagnostics 
    … …  
    866866/ 
    867867!----------------------------------------------------------------------- 
    868 &namtrd        !   diagnostics on dynamics and/or tracer trends         ("key_trddyn" and/or "key_trdtra") 
    869 !              !       or mixed-layer trends or barotropic vorticity    ("key_trdmld" or     "key_trdvor") 
    870 !----------------------------------------------------------------------- 
    871    nn_trd      = 365       !  time step frequency dynamics and tracers trends 
     868&namtrd        !   diagnostics on dynamics and/or tracer trends          
     869!              !       or mixed-layer trends or barotropic vorticity     
     870!----------------------------------------------------------------------- 
     871   ln_glo_trd  = .FALSE.   ! (T) global domain averaged diag for T, T^2, KE, and PE 
     872   ln_dyn_trd  = .FALSE.   ! (T) 3D momentum trend output 
     873   ln_dyn_mld  = .FALSE.   ! (T) 2D momentum trends averaged over the mixed layer  
     874   ln_vor_trd  = .FALSE.   ! (T) 2D barotropic vorticity trends  
     875   ln_KE_trd   = .FALSE.   ! (T) 3D Kinetic   Energy     trends 
     876   ln_PE_trd   = .FALSE.   ! (T) 3D Potential Energy     trends 
     877   ln_tra_trd  = .FALSE.   ! (T) 3D tracer trend output 
     878   ln_tra_mld  = .FALSE.   ! (T) 2D tracer trends averaged over the mixed layer  
     879   nn_trd      = 365       !  print frequency (ln_glo_trd=T) (unit=time step) 
     880/ 
     881!----------------------------------------------------------------------- 
     882&namtrd_mxl    !   mixed layer averaged trends diagnosed on dynamics and/or tracer trends          
     883!----------------------------------------------------------------------- 
     884   nn_trd      = 365       !  print frequency (ln_glo_trd=T) (unit=time step) 
     885   rn_zho_c    = 0.01      ! density criteria used to define the MLD (do not change unless changing the value used in zdfmxl) 
    872886   nn_ctls     =   0       !  control surface type in mixed-layer trends (0,1 or n<jpk) 
    873887   rn_ucf      =   1.      !  unit conversion factor (=1 -> /seconds ; =86400. -> /day) 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90

    r3294 r3876  
    1515   USE oce            ! ocean dynamics and tracers 
    1616   USE dom_oce        ! ocean space and time domain 
    17    USE trdmod_oce     ! ocean variables trends 
    18    USE trdmod         ! ocean dynamics trends 
     17   USE trd_oce        ! trends: ocean variables 
     18   USE trddyn         ! trend manager: dynamics 
    1919   USE in_out_manager ! I/O manager 
    2020   USE lib_mpp        ! MPP library 
    … …  
    103103         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 
    104104         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 
    105          CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt ) 
     105         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 
    106106         zfu_t(:,:,:) = ua(:,:,:) 
    107107         zfv_t(:,:,:) = va(:,:,:) 
    … …  
    153153         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 
    154154         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 
    155          CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt ) 
     155         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 
    156156      ENDIF 
    157157      !                                            ! Control print 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90

    r3294 r3876  
    1616   USE oce            ! ocean dynamics and tracers 
    1717   USE dom_oce        ! ocean space and time domain 
    18    USE trdmod         ! ocean dynamics trends 
    19    USE trdmod_oce     ! ocean variables trends 
     18   USE trd_oce        ! trends: ocean variables 
     19   USE trddyn         ! trend manager: dynamics 
    2020   USE in_out_manager ! I/O manager 
    2121   USE prtctl         ! Print control 
    … …  
    196196         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 
    197197         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 
    198          CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt ) 
     198         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 
    199199         zfu_t(:,:,:) = ua(:,:,:) 
    200200         zfv_t(:,:,:) = va(:,:,:) 
    … …  
    245245         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 
    246246         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 
    247          CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt ) 
     247         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 
    248248      ENDIF 
    249249      !                                            ! Control print 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r3294 r3876  
    1010 
    1111   !!---------------------------------------------------------------------- 
    12    !!   dyn_bfr      : Update the momentum trend with the bottom friction contribution 
     12   !!   dyn_bfr       : Update the momentum trend with the bottom friction contribution 
    1313   !!---------------------------------------------------------------------- 
    14    USE oce             ! ocean dynamics and tracers variables 
    15    USE dom_oce         ! ocean space and time domain variables  
    16    USE zdf_oce         ! ocean vertical physics variables 
    17    USE zdfbfr          ! ocean bottom friction variables 
    18    USE trdmod          ! ocean active dynamics and tracers trends  
    19    USE trdmod_oce      ! ocean variables trends 
    20    USE in_out_manager  ! I/O manager 
    21    USE prtctl          ! Print control 
    22    USE timing          ! Timing 
    23    USE wrk_nemo        ! Memory Allocation 
     14   USE oce            ! ocean dynamics and tracers variables 
     15   USE dom_oce        ! ocean space and time domain variables  
     16   USE zdf_oce        ! ocean vertical physics variables 
     17   USE zdfbfr         ! ocean bottom friction variables 
     18   USE trd_oce        ! trends: ocean variables 
     19   USE trddyn         ! trend manager: dynamics 
     20   USE in_out_manager ! I/O manager 
     21   USE prtctl         ! Print control 
     22   USE timing         ! Timing 
     23   USE wrk_nemo       ! Memory Allocation 
    2424 
    2525   IMPLICIT NONE 
    2626   PRIVATE 
    2727 
    28    PUBLIC   dyn_bfr    !  routine called by step.F90 
     28   PUBLIC   dyn_bfr   !  routine called by step.F90 
    2929 
    3030   !! * Substitutions 
    … …  
    5757      IF( nn_timing == 1 )  CALL timing_start('dyn_bfr') 
    5858      ! 
     59!!gm issue: better to put the logical in step to control the call of zdf_bfr 
     60!!          ==> change the logical from ln_bfrimp to ln_bfr_exp !! 
    5961      IF( .NOT.ln_bfrimp) THEN     ! only for explicit bottom friction form 
    6062                                    ! implicit bfr is implemented in dynzdf_imp 
    6163 
     64!!gm bug : time step is only rdt (not 2 rdt if euler start !) 
    6265        zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
    6366 
    … …  
    8992           ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    9093           ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    91            CALL trd_mod( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_trd_bfr, 'DYN', kt ) 
     94           CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 
    9295           CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    9396        ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r3764 r3876  
    3131   USE dom_oce         ! ocean space and time domain 
    3232   USE phycst          ! physical constants 
    33    USE trdmod          ! ocean dynamics trends 
    34    USE trdmod_oce      ! ocean variables trends 
     33   USE trd_oce         ! trends: ocean variables 
     34   USE trddyn          ! trend manager: dynamics 
    3535   USE in_out_manager  ! I/O manager 
    3636   USE prtctl          ! Print control 
    37    USE lbclnk          ! lateral boundary condition 
     37   USE lbclnk          ! lateral boundary condition  
    3838   USE lib_mpp         ! MPP library 
    3939   USE wrk_nemo        ! Memory Allocation 
    … …  
    4646   PUBLIC   dyn_hpg_init   ! routine called by opa module 
    4747 
    48    !                                              !!* Namelist namdyn_hpg : hydrostatic pressure gradient 
     48   !                                              !!* Namelist namdyn_hpg : hydrostatic pressure gradient  
    4949   LOGICAL , PUBLIC ::   ln_hpg_zco    = .TRUE.    !: z-coordinate - full steps 
    5050   LOGICAL , PUBLIC ::   ln_hpg_zps    = .FALSE.   !: z-coordinate - partial steps (interpolation) 
    … …  
    5454   LOGICAL , PUBLIC ::   ln_dynhpg_imp = .FALSE.   !: semi-implicite hpg flag 
    5555 
    56    INTEGER , PUBLIC ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
     56   INTEGER  ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 
    5757 
    5858   !! * Substitutions 
    … …  
    7070      !!                  ***  ROUTINE dyn_hpg  *** 
    7171      !! 
    72       !! ** Method  :   Call the hydrostatic pressure gradient routine 
     72      !! ** Method  :   Call the hydrostatic pressure gradient routine  
    7373      !!              using the scheme defined in the namelist 
    74       !! 
     74      !!    
    7575      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    76       !!             - Save the trend (l_trddyn=T) 
     76      !!             - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 
    7777      !!---------------------------------------------------------------------- 
    7878      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    … …  
    9999         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    100100         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    101          CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt ) 
     101         CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 
    102102         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    103103      ENDIF 
    … …  
    427427   END SUBROUTINE hpg_sco 
    428428 
     429 
    429430   SUBROUTINE hpg_djc( kt ) 
    430431      !!--------------------------------------------------------------------- 
    … …  
    664665      !! 
    665666      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    666       !!             - Save the trend (l_trddyn=T) 
    667       !! 
    668667      !!---------------------------------------------------------------------- 
    669668      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear 
    … …  
    717716 
    718717      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 
    719       DO jj = 1, jpj;   DO ji = 1, jpi 
    720           zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 
    721       END DO        ;   END DO 
    722  
    723       DO jk = 2, jpk;   DO jj = 1, jpj;   DO ji = 1, jpi 
    724           zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 
    725       END DO        ;   END DO        ;   END DO 
    726  
    727       fsp(:,:,:) = zrhh(:,:,:) 
     718      DO jj = 1, jpj 
     719         DO ji = 1, jpi 
     720            zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 
     721         END DO 
     722      END DO 
     723 
     724      DO jk = 2, jpk 
     725         DO jj = 1, jpj 
     726            DO ji = 1, jpi 
     727               zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 
     728            END DO 
     729         END DO 
     730      END DO 
     731 
     732      fsp(:,:,:) = zrhh (:,:,:) 
    728733      xsp(:,:,:) = zdept(:,:,:) 
    729734 
    … …  
    926931   END SUBROUTINE hpg_prj 
    927932 
     933 
    928934   SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 
    929935      !!---------------------------------------------------------------------- 
    … …  
    934940      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 
    935941      !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 
    936       !! 
    937942      !!---------------------------------------------------------------------- 
    938943      IMPLICIT NONE 
    … …  
    942947      INTEGER, INTENT(in) :: polynomial_type                        ! 1: cubic spline 
    943948                                                                    ! 2: Linear 
    944  
    945       ! Local Variables 
     949      ! 
    946950      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    947951      INTEGER  ::   jpi, jpj, jpkm1 
    … …  
    10331037      ENDIF 
    10341038 
    1035  
    10361039   END SUBROUTINE cspline 
    10371040 
    … …  
    10431046      !! ** Purpose :   1-d linear interpolation 
    10441047      !! 
    1045       !! ** Method  : 
    1046       !!                interpolation is straight forward 
     1048      !! ** Method  :   interpolation is straight forward 
    10471049      !!                extrapolation is also permitted (no value limit) 
    1048       !! 
    10491050      !!---------------------------------------------------------------------- 
    10501051      IMPLICIT NONE 
    … …  
    10631064   END FUNCTION interp1 
    10641065 
     1066 
    10651067   FUNCTION interp2(x, a, b, c, d)  RESULT(f) 
    10661068      !!---------------------------------------------------------------------- 
    … …  
    11261128   END FUNCTION integ_spline 
    11271129 
    1128  
    11291130   !!====================================================================== 
    11301131END MODULE dynhpg 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r3294 r3876  
    1414   USE oce             ! ocean dynamics and tracers 
    1515   USE dom_oce         ! ocean space and time domain 
    16    USE trdmod          ! ocean dynamics trends  
    17    USE trdmod_oce      ! ocean variables trends 
     16   USE trd_oce         ! trends: ocean variables 
     17   USE trddyn          ! trend manager: dynamics 
    1818   USE in_out_manager  ! I/O manager 
    1919   USE lib_mpp         ! MPP library 
    … …  
    5252      !! 
    5353      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 
    54       !!             - save this trends (l_trddyn=T) for post-processing 
     54      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing 
    5555      !!---------------------------------------------------------------------- 
    5656      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    … …  
    131131         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    132132         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    133          CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt ) 
     133         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 
    134134         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    135135      ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

    r3294 r3876  
    2020   USE dynldf_iso     ! lateral mixing            (dyn_ldf_iso    routine) 
    2121   USE dynldf_lap     ! lateral mixing            (dyn_ldf_lap    routine) 
    22    USE trdmod         ! ocean dynamics and tracer trends 
    23    USE trdmod_oce     ! ocean variables trends 
     22   USE trd_oce        ! trends: ocean variables 
     23   USE trddyn         ! trend manager: dynamics 
    2424   USE prtctl         ! Print control 
    2525   USE in_out_manager ! I/O manager 
    … …  
    5454      !! ** Purpose :   compute the lateral ocean dynamics physics. 
    5555      !!---------------------------------------------------------------------- 
    56       ! 
    5756      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    5857      ! 
    … …  
    106105         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    107106         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    108          CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_ldf, 'DYN', kt ) 
     107         CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 
    109108         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
    110109      ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90

    r3634 r3876  
    1919   USE dom_oce         ! ocean space and time domain 
    2020   USE ldfdyn_oce      ! ocean dynamics: lateral physics 
     21   ! 
    2122   USE in_out_manager  ! I/O manager 
    22    USE trdmod          ! ocean dynamics trends  
    23    USE trdmod_oce      ! ocean variables trends 
    2423   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2524   USE wrk_nemo        ! Memory Allocation 
    … …  
    7069      !!      Add this before trend to the general trend (ua,va): 
    7170      !!            (ua,va) = (ua,va) + (diffu,diffv) 
    72       !!      'key_trddyn' defined: the two components of the horizontal 
    73       !!                               diffusion trend are saved. 
    7471      !! 
    7572      !! ** Action : - Update (ua,va) with the before iso-level biharmonic 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90

    r3634 r3876  
    2020   USE ldfdyn_oce      ! ocean dynamics lateral physics 
    2121   USE zdf_oce         ! ocean vertical physics 
    22    USE trdmod          ! ocean dynamics trends  
    23    USE trdmod_oce      ! ocean variables trends 
    2422   USE ldfslp          ! iso-neutral slopes available 
     23   ! 
    2524   USE in_out_manager  ! I/O manager 
    2625   USE lib_mpp         ! MPP library 
    … …  
    8180      !!         -3- Add this trend to the general trend (ta,sa): 
    8281      !!            (ua,va) = (ua,va) + (zwk3,zwk4) 
    83       !!      'key_trddyn' defined: the trend is saved for diagnostics. 
    8482      !! 
    8583      !! ** Action  : - Update (ua,va) arrays with the before geopotential 
    8684      !!                biharmonic mixing trend. 
    87       !!              - save the trend in (zwk3,zwk4) ('key_trddyn') 
    8885      !!---------------------------------------------------------------------- 
    8986      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index 
    … …  
    176173      !!                          pu and pv (all the components except 
    177174      !!                          second order vertical derivative term) 
    178       !!      'key_trddyn' defined: the trend is saved for diagnostics. 
    179175      !!---------------------------------------------------------------------- 
    180176      !! 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

    r3294 r3876  
    2222   USE ldftra_oce      ! ocean tracer   lateral physics 
    2323   USE zdf_oce         ! ocean vertical physics 
    24    USE trdmod          ! ocean dynamics trends  
    25    USE trdmod_oce      ! ocean variables trends 
    2624   USE ldfslp          ! iso-neutral slopes  
    2725   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90

    r3294 r3876  
    1919   USE ldfdyn_oce      ! ocean dynamics: lateral physics 
    2020   USE zdf_oce         ! ocean vertical physics 
     21!!gm ??   USE ldfslp          ! iso-neutral slopes  
     22   ! 
    2123   USE in_out_manager  ! I/O manager 
    22    USE trdmod          ! ocean dynamics trends  
    23    USE trdmod_oce      ! ocean variables trends 
    24    USE ldfslp          ! iso-neutral slopes  
    2524   USE timing          ! Timing 
    2625 
    … …  
    5756      !!      Add this before trend to the general trend (ua,va): 
    5857      !!            (ua,va) = (ua,va) + (diffu,diffv) 
    59       !!      'key_trddyn' activated: the two components of the horizontal 
    60       !!                                 diffusion trend are saved. 
    6158      !! 
    6259      !! ** Action : - Update (ua,va) with the before iso-level harmonic  

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90

    r3723 r3876  
    55   !!                 recoded version of simplest case (u*, v* only) 
    66   !!====================================================================== 
    7    !! History :  1.0  !  2007-06  (Michael Dunphy)  Modular form: - new namelist parameters 
    8    !!                                                             - horizontal diffusion for Neptune 
    9    !!                                                             - vertical diffusion for gm in momentum eqns 
    10    !!                                                             - option to use Neptune in Coriolis eqn 
     7   !! History :  1.0  !  2007-06  (Zeliang Wang, Michael Dunphy, BIO)  Original code 
     8   !!                             Modular form: - new namelist parameters 
     9   !!                                           - horizontal diffusion for Neptune 
     10   !!                                           - vertical diffusion for gm in momentum eqns 
     11   !!                                           - option to use Neptune in Coriolis eqn 
    1112   !!                    2011-08  (Jeff Blundell, NOCS) Simplified form for temporally invariant u*, v* 
    12    !!                                               Horizontal and vertical diffusivity formulations removed 
    13    !!                                               Dynamic allocation of storage added 
    14    !!                                               Option of ramping Neptune vel. down 
    15    !!                                               to zero added in shallow depths added 
     13   !!                             Horizontal and vertical diffusivity formulations removed 
     14   !!                             Dynamic allocation of storage added 
     15   !!                             Option of ramping Neptune vel. down to zero added in shallow depths added 
     16   !!---------------------------------------------------------------------- 
     17 
    1618   !!---------------------------------------------------------------------- 
    1719   !! dynnept_alloc        : 
    … …  
    3032   USE phycst 
    3133   USE lbclnk 
    32    USE wrk_nemo        ! Memory Allocation 
     34   USE wrk_nemo         ! Memory Allocation 
    3335 
    3436   IMPLICIT NONE 
    3537   PRIVATE 
    3638 
    37    !! * Routine accessibility 
    38    PUBLIC dyn_nept_init      ! routine called by nemogcm.F90 
    39    PUBLIC dyn_nept_cor       ! routine called by step.F90 
    40    !! dynnept_alloc()         is called only by dyn_nept_init, within this module 
    41    !! dyn_nept_div_cur_init   is called only by dyn_nept_init, within this module 
    42    !! dyn_nept_vel            is called only by dyn_nept_cor,  within this module 
    43  
    44    !! * Shared module variables 
    45    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: zunep, zvnep  ! Neptune u and v 
    46    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: zhdivnep      ! hor. div for Neptune vel. 
    47    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: zmrotnep      ! curl for Neptune vel. 
    48  
    49  
    50    !! * Namelist namdyn_nept variables 
     39   PUBLIC   dyn_nept_init   ! routine called by nemogcm.F90 
     40   PUBLIC   dyn_nept_cor    ! routine called by step.F90 
     41 
     42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zunep, zvnep   ! Neptune u and v 
     43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zhdivnep       ! hor. div for Neptune vel. 
     44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zmrotnep       ! curl for Neptune vel. 
     45 
     46   !                                                 !!* Namelist namdyn_nept variables 
    5147   LOGICAL, PUBLIC  ::  ln_neptsimp        = .FALSE.  ! yes/no simplified neptune 
    5248 
    … …  
    6056   REAL(wp)         ::  rn_htrmax          =  200.0   ! max. depth of transition range 
    6157 
    62    !! * Module variables 
    63  
    64  
    6558   !! * Substitutions 
    6659#  include "vectopt_loop_substitute.h90" 
    6760#  include "domzgr_substitute.h90" 
    6861   !!---------------------------------------------------------------------- 
    69    !!   OPA 9.0 , implemented by Bedford Institute of Oceanography 
    70    !!---------------------------------------------------------------------- 
    71  
     62   !! NEMO/OPA 3.3 , NEMO Consortium (2011) 
     63   !! $Id: dynadv_cen2.F90 3316 2012-02-21 16:00:02Z gm $ 
     64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     65   !!---------------------------------------------------------------------- 
    7266 CONTAINS 
    7367 
    … …  
    9084      !!                and compute the arrays zunep and zvnep 
    9185      !! 
    92       !! ** Method  :   zunep = 
    93       !!                zvnep = 
    94       !! 
    95       !! ** History :  1.0  !   07-05  (Zeliang Wang)   Original code for zunep, zvnep 
    96       !!               1.1  !   07-06  (Michael Dunphy) namelist and  initialisation 
    97       !!               2.0  ! 2011-07  (Jeff Blundell, NOCS) 
    98       !!                    ! Simplified form for temporally invariant u*, v* 
    99       !!                    ! Horizontal and vertical diffusivity formulations removed 
    100       !!                    ! Includes optional tapering-off in shallow depths 
     86      !! ** Method  :   Simplified form for temporally invariant u*, v* 
     87      !!                Horizontal and vertical diffusivity formulations removed 
     88      !!                Includes optional tapering-off in shallow depths 
    10189      !!---------------------------------------------------------------------- 
    10290      USE iom 
    103       !! 
     91      ! 
    10492      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    10593      REAL(wp) :: unemin,unemax,vnemin,vnemax   ! extrema of (u*, v*) fields 
    … …  
    140128      ENDIF 
    141129      ! 
    142       IF( .NOT. ln_neptsimp ) RETURN 
     130      IF( .NOT. ln_neptsimp )   RETURN 
    143131      !                                 ! Dynamically allocate local work arrays 
    144132      CALL wrk_alloc( jpi, jpj     , ht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n  )  
    … …  
    352340      !! ** Action  : - compute zhdivnep, the hor. divergence of (u*, v*) 
    353341      !!              - compute zmrotnep, the rel. vorticity  of (u*, v*) 
    354       !! 
    355       !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code 
    356       !!            4.0  ! 1991-11  (G. Madec) 
    357       !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions 
    358       !!            7.0  ! 1996-01  (G. Madec)  s-coordinates 
    359       !!            8.0  ! 1997-06  (G. Madec)  lateral boundary cond., lbc 
    360       !!            8.1  ! 1997-08  (J.M. Molines)  Open boundaries 
    361       !!            8.2  ! 2000-03  (G. Madec)  no slip accurate 
    362       !!  NEMO      1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90 
    363       !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries 
    364       !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90 
    365       !!             -   ! 2005-01  (J. Chanut, A. Sellar) unstructured open boundaries 
    366       !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    367       !!                 ! 2011-06  (Jeff Blundell, NOCS) Adapt code from divcur.F90 
    368       !!                 !           to compute Neptune effect fields only 
    369       !!---------------------------------------------------------------------- 
    370       ! 
     342      !!---------------------------------------------------------------------- 
    371343      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
    372344      !!---------------------------------------------------------------------- 
    … …  
    508480 
    509481   SUBROUTINE dyn_nept_smooth_vel( htold, htnew, ld_option ) 
    510  
    511482      !!---------------------------------------------------------------------- 
    512483      !!                  ***  ROUTINE dyn_nept_smooth_vel  *** 
    … …  
    598569   END SUBROUTINE dyn_nept_smooth_vel 
    599570 
     571   !!============================================================================== 
    600572END MODULE dynnept 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r3764 r3876  
    1717   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module 
    1818   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 
     19   !!            3.4  !  2012-02  (G. Madec) add the diagnostic of the time filter trends 
    1920   !!------------------------------------------------------------------------- 
    2021   
    … …  
    2930   USE dynadv          ! dynamics: vector invariant versus flux form 
    3031   USE domvvl          ! variable volume 
     32   USE trd_oce         ! trends: ocean variables 
     33   USE trddyn          ! trend manager: dynamics 
     34   USE trdken          ! trend manager: kinetic energy 
    3135   USE obc_oce         ! ocean open boundary conditions 
    3236   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine) 
    … …  
    3741   USE bdydyn          ! ocean open boundary conditions 
    3842   USE bdyvol          ! ocean open boundary condition (bdy_vol routines) 
     43   ! 
    3944   USE in_out_manager  ! I/O manager 
     45   USE iom             ! I/O manager library 
    4046   USE lbclnk          ! lateral boundary condition (or mpp link) 
    4147   USE lib_mpp         ! MPP library 
    4248   USE wrk_nemo        ! Memory Allocation 
    4349   USE prtctl          ! Print control 
     50   USE timing          ! Timing 
    4451#if defined key_agrif 
    4552   USE agrif_opa_interp 
    4653#endif 
    47    USE timing          ! Timing 
    4854 
    4955   IMPLICIT NONE 
    … …  
    8187      !!             at the local domain boundaries through lbc_lnk call, 
    8288      !!             at the one-way open boundaries (lk_obc=T), 
    83       !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
     89      !!             at the AGRIF zoom   boundaries (lk_agrif=T) 
    8490      !! 
    8591      !!              * Apply the time filter applied and swap of the dynamics 
    … …  
    101107      REAL(wp) ::   z2dt         ! temporary scalar 
    102108#endif 
    103       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars 
    104       REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      - 
    105       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f  
     109      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars 
     110      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
     111      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva  
    106112      !!---------------------------------------------------------------------- 
    107113      ! 
    108       IF( nn_timing == 1 )  CALL timing_start('dyn_nxt') 
    109       ! 
    110       CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     114      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt') 
     115      ! 
     116      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 
    111117      ! 
    112118      IF( kt == nit000 ) THEN 
    … …  
    156162# if defined key_obc 
    157163      !                                !* OBC open boundaries 
    158       IF( lk_obc ) CALL obc_dyn( kt ) 
     164      IF( lk_obc )   CALL obc_dyn( kt ) 
    159165      ! 
    160166      IF( .NOT. lk_dynspg_flt ) THEN 
    … …  
    163169         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T) 
    164170         !                              - Correct the barotropic velocities 
    165          IF( lk_obc ) CALL obc_dyn_bt( kt ) 
     171         IF( lk_obc )   CALL obc_dyn_bt( kt ) 
    166172         ! 
    167173!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called 
    … …  
    186192# endif 
    187193#endif 
     194 
     195    IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics 
     196         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step  
     197         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt 
     198         ! 
     199         !                                  ! Kinetic energy and Conversion 
     200         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt ) 
     201         ! 
     202         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends 
     203            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 
     204            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 
     205            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter 
     206            CALL iom_put( "vtrd_tot", zva ) 
     207         ENDIF 
     208         ! 
     209         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter 
     210         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the 
     211         !                                  !  computation of the asselin filter trends) 
     212      ENDIF 
    188213 
    189214      ! Time filter and swap of dynamics arrays 
    … …  
    274299      ENDIF 
    275300 
     301      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum 
     302         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 
     303         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 
     304         CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 
     305      ENDIF 
     306 
    276307      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    277308         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    278309      !  
    279       CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     310      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 
    280311      ! 
    281312      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt') 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r3625 r3876  
    2222   USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine) 
    2323   USE dynadv         ! dynamics: vector invariant versus flux form 
    24    USE trdmod         ! ocean dynamics trends 
    25    USE trdmod_oce     ! ocean variables trends 
     24   USE trd_oce        ! trends: ocean variables 
     25   USE trddyn         ! trend manager: dynamics 
    2626   USE prtctl         ! Print control                     (prt_ctl routine) 
    2727   USE in_out_manager ! I/O manager 
    2828   USE lib_mpp        ! MPP library 
    29    USE solver          ! solver initialization 
    30    USE wrk_nemo        ! Memory Allocation 
    31    USE timing          ! Timing 
     29   USE solver         ! solver initialization 
     30   USE wrk_nemo       ! Memory Allocation 
     31   USE timing         ! Timing 
    3232 
    3333 
    … …  
    118118            END DO 
    119119         END DO 
     120          
     121!!gm add here a call to dyn_trd for apr, the surf pressure trends ???? 
     122          
    120123      ENDIF 
    121124 
    … …  
    142145         ! 
    143146         CALL wrk_dealloc( jpi, jpj, zpice ) 
     147          
     148!!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ???? 
     149          
    144150      ENDIF 
    145151 
    … …  
    170176         CASE( 2 ) 
    171177            z2dt = 2. * rdt 
    172             IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 
     178            IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
    173179            ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 
    174180            ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 
    175181         END SELECT 
    176          CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_spg, 'DYN', kt ) 
     182         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 
    177183         ! 
    178184         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

    r3680 r3876  
    2222   USE obc_par         ! open boundary condition parameters 
    2323   USE obcdta          ! open boundary condition data     (bdy_dta_bt routine) 
     24   ! 
    2425   USE in_out_manager  ! I/O manager 
    2526   USE lib_mpp         ! distributed memory computing library 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r3765 r3876  
    1313   !!             -   !  2006-08  (J.Chanut, A.Sellar) Calls to BDY routines.  
    1414   !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 
     15   !!            3.5  !  2012-05  (F. Roquet, G. Madec)  add some trends diag 
    1516   !!---------------------------------------------------------------------- 
    1617#if defined key_dynspg_flt   ||   defined key_esopa   
    … …  
    3940   USE bdyvol          ! ocean open boundary condition (bdy_vol routine) 
    4041   USE cla             ! cross land advection 
     42   USE trd_oce         ! trends: ocean variables 
     43   USE trddyn          ! trend manager: dynamics 
     44   ! 
    4145   USE in_out_manager  ! I/O manager 
    4246   USE lib_mpp         ! distributed memory computing library 
    … …  
    4650   USE iom 
    4751   USE lib_fortran 
     52   USE timing          ! Timing 
    4853#if defined key_agrif 
    4954   USE agrif_opa_interp 
    5055#endif 
    51    USE timing          ! Timing 
    5256 
    5357   IMPLICIT NONE 
    … …  
    102106      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 
    103107      !! 
    104       !! References : Roullet and Madec 1999, JGR. 
     108      !! References : Roullet and Madec, JGR, 2000. 
    105109      !!--------------------------------------------------------------------- 
    106110      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    107111      INTEGER, INTENT(  out) ::   kindic   ! solver convergence flag (<0 if not converge) 
    108       !!                                    
     112      ! 
    109113      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    110114      REAL(wp) ::   z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv   ! local scalars 
    111       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zub, zvb 
     115      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     116      REAL(wp), POINTER, DIMENSION(:,:)   ::  zpw 
    112117      !!---------------------------------------------------------------------- 
    113118      ! 
    114119      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt') 
    115       ! 
    116       CALL wrk_alloc( jpi,jpj,jpk, zub, zvb ) 
    117120      ! 
    118121      IF( kt == nit000 ) THEN 
    … …  
    183186            END DO 
    184187         END DO 
     188         ! save the explicit SPG trends for further diagnostics 
     189         ! --------------------------------------------- 
     190         IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
     191            CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
     192            DO jk = 1, jpkm1              ! unweighted time stepping  
     193               DO jj = 2, jpjm1 
     194                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     195                     ztrdu(ji,jj,jk) = spgu(ji,jj) * umask(ji,jj,jk) 
     196                     ztrdv(ji,jj,jk) = spgv(ji,jj) * vmask(ji,jj,jk) 
     197                  END DO 
     198               END DO 
     199            END DO 
     200            CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgexp, kt ) 
     201         ENDIF 
    185202         ! 
    186203      ENDIF 
    … …  
    208225      END DO 
    209226 
    210       ! vertical sum 
     227      !                                    ! vertical sum 
    211228!CDIR NOLOOPCHG 
    212229      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    … …  
    217234            END DO 
    218235         END DO 
    219       ELSE                        ! No  vector opt. 
     236      ELSE                                      ! No  vector opt. 
    220237         DO jk = 1, jpkm1 
    221238            DO jj = 2, jpjm1 
    … …  
    228245      ENDIF 
    229246 
    230       ! transport: multiplied by the horizontal scale factor 
    231       DO jj = 2, jpjm1 
     247      DO jj = 2, jpjm1                     ! transport: multiplied by the horizontal scale factor 
    232248         DO ji = fs_2, fs_jpim1   ! vector opt. 
    233249            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj) 
    … …  
    342358      ENDIF 
    343359#endif       
     360 
     361      IF( l_trddyn )   THEN                      
     362         ztrdu(:,:,:) = ua(:,:,:)                 ! save the after velocity before the filtered SPG 
     363         ztrdv(:,:,:) = va(:,:,:) 
     364         ! 
     365         CALL wrk_alloc( jpi, jpj, zpw ) 
     366         ! 
     367         zpw(:,:) = - z2dt * gcx(:,:) 
     368         CALL iom_put( "ssh_flt" , zpw )          ! output equivalent ssh modification due to implicit filter 
     369         ! 
     370         !                                        ! save surface pressure flux: -pw at z=0 
     371         zpw(:,:) = - rau0 * grav * sshn(:,:) * wn(:,:,1) * tmask(:,:,1) 
     372         CALL iom_put( "pw0_exp" , zpw ) 
     373         zpw(:,:) = wn(:,:,1) 
     374         CALL iom_put( "w0" , zpw ) 
     375         zpw(:,:) =  rau0 * z2dtg * gcx(:,:) * wn(:,:,1) * tmask(:,:,1) 
     376         CALL iom_put( "pw0_flt" , zpw ) 
     377         ! 
     378         CALL wrk_dealloc( jpi, jpj, zpw )  
     379         !                                    
     380      ENDIF 
     381       
    344382      ! Add the trends multiplied by z2dt to the after velocity 
    345383      ! ------------------------------------------------------- 
    … …  
    356394      END DO 
    357395 
    358       ! write filtered free surface arrays in restart file 
    359       ! -------------------------------------------------- 
    360       IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 
    361       ! 
    362       CALL wrk_dealloc( jpi,jpj,jpk, zub, zvb ) 
     396      IF( l_trddyn )   THEN                      ! save the explicit SPG trends for further diagnostics 
     397         ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / z2dt 
     398         ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / z2dt 
     399         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgflt, kt ) 
     400         ! 
     401         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  
     402      ENDIF 
     403 
     404      IF( lrst_oce )   CALL flt_rst( kt, 'WRITE' )      ! write filtered free surface arrays in restart file 
    363405      ! 
    364406      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt') 
    … …  
    373415      !! ** Purpose : Read or write filtered free surface arrays in restart file 
    374416      !!---------------------------------------------------------------------- 
    375       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    376       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
     417      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     418      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    377419      !!---------------------------------------------------------------------- 
    378420      ! 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r3680 r3876  
    400400         !                                                !* Update the forcing (BDY and tides) 
    401401         !                                                !  ------------------ 
    402          IF( lk_obc )   CALL obc_dta_bt ( kt, jn   ) 
    403          IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 
    404          IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn ) 
     402         IF( lk_obc                    )   CALL obc_dta_bt ( kt, jn   ) 
     403         IF( lk_bdy                    )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 
     404         IF( ln_tide_pot .AND. lk_tide )  CALL upd_tide( kt, jn ) 
    405405 
    406406         !                                                !* after ssh_e 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r3802 r3876  
    1515   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme 
    1616   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     17   !!            3.5  ! 2012-02  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity  
    1718   !!---------------------------------------------------------------------- 
    1819 
    … …  
    2930   USE dommsk         ! ocean mask 
    3031   USE dynadv         ! momentum advection (use ln_dynadv_vec value) 
    31    USE trdmod         ! ocean dynamics trends  
    32    USE trdmod_oce     ! ocean variables trends 
     32   USE trd_oce        ! trends: ocean variables 
     33   USE trddyn         ! trend manager: dynamics 
    3334   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3435   USE prtctl         ! Print control 
    … …  
    7374      !! ** Action : - Update (ua,va) with the now vorticity term trend 
    7475      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    75       !!               and planetary vorticity trends) ('key_trddyn') 
     76      !!               and planetary vorticity trends) (l_trddyn=T) 
    7677      !!---------------------------------------------------------------------- 
    7778      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    … …  
    108109            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    109110            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    110             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt ) 
     111            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
    111112            ztrdu(:,:,:) = ua(:,:,:) 
    112113            ztrdv(:,:,:) = va(:,:,:) 
    … …  
    114115            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    115116            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    116             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 
    117             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 
     117            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    118118         ELSE 
    119119            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity 
    … …  
    127127            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    128128            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    129             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt ) 
     129            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
    130130            ztrdu(:,:,:) = ua(:,:,:) 
    131131            ztrdv(:,:,:) = va(:,:,:) 
    … …  
    133133            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    134134            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    135             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 
    136             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 
     135            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    137136         ELSE 
    138137            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity 
    … …  
    146145            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    147146            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    148             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt ) 
     147            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
    149148            ztrdu(:,:,:) = ua(:,:,:) 
    150149            ztrdv(:,:,:) = va(:,:,:) 
    … …  
    152151            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    153152            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    154             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 
    155             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 
     153            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    156154         ELSE 
    157155            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene) 
    … …  
    165163            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    166164            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    167             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt ) 
     165            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
    168166            ztrdu(:,:,:) = ua(:,:,:) 
    169167            ztrdv(:,:,:) = va(:,:,:) 
    … …  
    171169            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    172170            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    173             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 
    174             CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 
     171            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    175172         ELSE 
    176173            CALL vor_een( kt, ntot, ua, va )                ! total vorticity 
    … …  
    211208      !! 
    212209      !! ** Action : - Update (ua,va) with the now vorticity term trend 
    213       !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    214       !!               and planetary vorticity trends) ('key_trddyn') 
    215210      !! 
    216211      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    … …  
    328323      !! 
    329324      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 
    330       !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    331       !!               and planetary vorticity trends) ('key_trddyn') 
    332325      !! 
    333326      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    … …  
    444437      !! 
    445438      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 
    446       !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative  
    447       !!               and planetary vorticity trends) ('key_trddyn') 
    448439      !! 
    449440      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    … …  
    557548      !! 
    558549      !! ** Action : - Update (ua,va) with the now vorticity term trend 
    559       !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    560       !!               and planetary vorticity trends) ('key_trddyn') 
    561550      !! 
    562551      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r3294 r3876  
    1616   USE dom_oce        ! ocean space and time domain 
    1717   USE sbc_oce        ! surface boundary condition: ocean 
    18    USE trdmod_oce     ! ocean variables trends 
    19    USE trdmod         ! ocean dynamics trends  
     18   USE trd_oce        ! trends: ocean variables 
     19   USE trddyn         ! trend manager: dynamics 
     20   ! 
    2021   USE in_out_manager ! I/O manager 
    21    USE lib_mpp         ! MPP library 
     22   USE lib_mpp        ! MPP library 
    2223   USE prtctl         ! Print control 
    23    USE wrk_nemo        ! Memory Allocation 
    24    USE timing          ! Timing 
     24   USE wrk_nemo       ! Memory Allocation 
     25   USE timing         ! Timing 
    2526 
    2627   IMPLICIT NONE 
    … …  
    5354      !! 
    5455      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends 
    55       !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn') 
     56      !!              - Send the trends to trddyn for diagnostics (l_trddyn=T) 
    5657      !!---------------------------------------------------------------------- 
    5758      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx 
    … …  
    118119         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    119120         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    120          CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt) 
     121         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 
    121122         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    122123      ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

    r3294 r3876  
    2020 
    2121   USE ldfdyn_oce      ! ocean dynamics: lateral physics 
    22    USE trdmod          ! ocean active dynamics and tracers trends  
    23    USE trdmod_oce      ! ocean variables trends 
     22   USE trd_oce         ! trends: ocean variables 
     23   USE trddyn          ! trend manager: dynamics 
    2424   USE in_out_manager  ! I/O manager 
    2525   USE lib_mpp         ! MPP library 
    … …  
    9191         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    9292         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    93          CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_zdf, 'DYN', kt ) 
     93         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 
    9494         ! 
    9595         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r3625 r3876  
    9494      IF( ln_bfrimp ) THEN 
    9595# if defined key_vectopt_loop 
    96       DO jj = 1, 1 
    97          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     96         DO jj = 1, 1 
     97            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    9898# else 
    99       DO jj = 2, jpjm1 
    100          DO ji = 2, jpim1 
     99         DO jj = 2, jpjm1 
     100            DO ji = 2, jpim1 
    101101# endif 
    102             ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    103             ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    104             zavmu(ji,jj) = avmu(ji,jj,ikbu+1) 
    105             zavmv(ji,jj) = avmv(ji,jj,ikbv+1) 
    106             avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)  
    107             avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
    108          END DO 
    109       END DO 
     102               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     103               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     104               zavmu(ji,jj) = avmu(ji,jj,ikbu+1) 
     105               zavmv(ji,jj) = avmv(ji,jj,ikbv+1) 
     106               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)  
     107               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
     108            END DO 
     109         END DO 
    110110      ENDIF 
    111111 
    … …  
    284284      IF( ln_bfrimp ) THEN 
    285285# if defined key_vectopt_loop 
    286       DO jj = 1, 1 
    287          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     286         DO jj = 1, 1 
     287            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    288288# else 
    289       DO jj = 2, jpjm1 
    290          DO ji = 2, jpim1 
     289         DO jj = 2, jpjm1 
     290            DO ji = 2, jpim1 
    291291# endif 
    292             ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    293             ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    294             avmu(ji,jj,ikbu+1) = zavmu(ji,jj) 
    295             avmv(ji,jj,ikbv+1) = zavmv(ji,jj) 
    296          END DO 
    297       END DO 
     292               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     293               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     294               avmu(ji,jj,ikbu+1) = zavmu(ji,jj) 
     295               avmv(ji,jj,ikbv+1) = zavmv(ji,jj) 
     296            END DO 
     297         END DO 
    298298      ENDIF 
    299299      ! 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r3625 r3876  
    1818   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
    1919   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp 
     20   !!            3.5  ! 2012-03  (F. Roquet, G. Madec)  add pen_ddt_dds used in trdpen 
     21   !!             -   ! 2012-05  (F. Roquet)  add Vallis and original JM95 equation of state 
     22   !!             -   ! 2012-05  (F. Roquet)  add eos_alpha_beta 
    2023   !!---------------------------------------------------------------------- 
    2124 
    … …  
    2831   !!   eos_bn2        : Compute the Brunt-Vaisala frequency 
    2932   !!   eos_alpbet     : calculates the in situ thermal/haline expansion ratio 
     33   !!   eos_expansion_ratio : calculates the partial derivative of in situ density with respect to T and S 
    3034   !!   tfreez         : Compute the surface freezing temperature 
    3135   !!   eos_init       : set eos parameters (namelist) 
    … …  
    5155   END INTERFACE 
    5256 
    53    PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules 
    54    PUBLIC   eos_init   ! called by istate module 
    55    PUBLIC   bn2        ! called by step module 
    56    PUBLIC   eos_alpbet ! called by ldfslp module 
    57    PUBLIC   tfreez     ! called by sbcice_... modules 
     57   PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules 
     58   PUBLIC   eos_init       ! called by istate module 
     59   PUBLIC   bn2            ! called by step module 
     60   PUBLIC   eos_alpbet     ! called by ldfslp module 
     61   PUBLIC   eos_alpha_beta ! used for density diagnostics in dyntra 
     62   PUBLIC   eos_pen        ! used for pe diagnostics in trdpen 
     63   PUBLIC   tfreez         ! called by sbcice_... modules 
    5864 
    5965   !                                          !!* Namelist (nameos) * 
    60    INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
     66   INTEGER , PUBLIC ::   nn_eos   = 0         !: = -1/0/1/2/3 type of eq. of state and associated Brunt-Vaisala frequ. 
    6167   REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state) 
    6268   REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4_wp !: saline  expension coeff. (linear equation of state) 
    6369 
    64    REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio 
     70   REAL(wp) ::   vallis_ratio = 0   ! 1027/rau0 
     71   REAL(wp) ::   vallis_diff  = 0   ! (1027-rau0)/rau0 
    6572 
    6673   !! * Substitutions 
    … …  
    8289      !!       defined through the namelist parameter nn_eos. 
    8390      !! 
    84       !! ** Method  :   3 cases: 
    85       !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state. 
     91      !! ** Method  :   5 cases: 
     92      !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
     93      !!         Check value: rho = 1041.83267 kg/m**3 for p=3000 dbar, 
     94      !!          t = 3 deg celcius, s=35.5 psu 
     95      !!      nn_eos = 0 : standard NEMO equation of state, based on  
     96      !!         Jackett and McDougall (1995) equation of state. 
    8697      !!         the in situ density is computed directly as a function of 
    8798      !!         potential temperature relative to the surface (the opa t 
    … …  
    103114      !!               salinity 
    104115      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 
     116      !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 
     117      !!              prd(t,s,p) =  ( rho(t,s,p) - rau0 ) / rau0 
     118      !!              where the pressure p in decibars is approximated by the depth in meters. 
    105119      !!      Note that no boundary condition problem occurs in this routine 
    106120      !!      as pts are defined over the whole domain. 
    … …  
    108122      !! ** Action  :   compute prd , the in situ density (no units) 
    109123      !! 
    110       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 
    111       !!---------------------------------------------------------------------- 
    112       !! 
     124      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     125      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
     126      !!---------------------------------------------------------------------- 
    113127      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    114128      !                                                      ! 2 : salinity               [psu] 
    115129      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-] 
    116       !! 
     130      ! 
    117131      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    118132      REAL(wp) ::   zt , zs , zh , zsr   ! local scalars 
    … …  
    123137      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    124138      !!---------------------------------------------------------------------- 
    125  
    126139      ! 
    127140      IF( nn_timing == 1 ) CALL timing_start('eos') 
    … …  
    131144      SELECT CASE( nn_eos ) 
    132145      ! 
    133       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
     146      CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==! 
     147         ! 
    134148!CDIR NOVERRCHK 
    135149         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    … …  
    156170                  ! 
    157171                  ! add the compression terms 
     172                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     173                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     174                  zb = zbw + ze * zs 
     175                  ! 
     176                  zd = -1.480266e-4_wp 
     177                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
     178                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
     179                  za = ( zd*zsr + zc ) *zs + zaw 
     180                  ! 
     181                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     182                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     183                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
     184                  zk0= ( zb1*zsr + za1 )*zs + zkw 
     185                  ! 
     186                  ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 
     187                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
     188                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
     189               END DO 
     190            END DO 
     191         END DO 
     192         ! 
     193      CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==! 
     194         ! 
     195!CDIR NOVERRCHK 
     196         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     197         ! 
     198         DO jk = 1, jpkm1 
     199            DO jj = 1, jpj 
     200               DO ji = 1, jpi 
     201                  zt = pts   (ji,jj,jk,jp_tem) 
     202                  zs = pts   (ji,jj,jk,jp_sal) 
     203                  zh = fsdept(ji,jj,jk)        ! depth 
     204                  zsr= zws   (ji,jj,jk)        ! square root salinity 
     205                  ! 
     206                  ! compute volumic mass pure water at atm pressure 
     207                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     208                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
     209                  ! seawater volumic mass atm pressure 
     210                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     211                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     212                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     213                  zr4= 4.8314e-4_wp 
     214                  ! 
     215                  ! potential volumic mass (reference to the surface) 
     216                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     217                  ! 
     218                  ! add the compression terms 
    158219                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    159220                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    … …  
    187248         END DO 
    188249         ! 
     250      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     251         DO jk = 1, jpkm1 
     252            DO jj = 1, jpj 
     253               DO ji = 1, jpi 
     254                  zt = pts   (ji,jj,jk,jp_tem) - 10._wp 
     255                  zs = pts   (ji,jj,jk,jp_sal) - 35._wp 
     256                  zh = fsdept(ji,jj,jk)        ! depth 
     257                  ! masked in situ density anomaly 
     258                  prd(ji,jj,jk) = vallis_diff + vallis_ratio * (                                             & 
     259                     &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  & 
     260                     &   ) * tmask(ji,jj,jk) 
     261               END DO 
     262            END DO 
     263         END DO 
     264         ! 
    189265      END SELECT 
    190266      ! 
    … …  
    207283      !!     namelist parameter nn_eos. 
    208284      !! 
    209       !! ** Method  : 
    210       !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state. 
     285      !! ** Method  :   5 cases: 
     286      !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
     287      !!      nn_eos = 0 : standard NEMO equation of state, based on  
     288      !!         Jackett and McDougall (1995) equation of state. 
    211289      !!         the in situ density is computed directly as a function of 
    212290      !!         potential temperature relative to the surface (the opa t 
    … …  
    235313      !!                       = rn_beta * s - rn_alpha * tn - 1. 
    236314      !!              rhop(t,s)  = rho(t,s) 
     315      !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 
     316      !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 
     317      !!              rhop(t,s)  = rho(t,s,0) 
    237318      !!      Note that no boundary condition problem occurs in this routine 
    238319      !!      as (tn,sn) or (ta,sa) are defined over the whole domain. 
    … …  
    241322      !!              - prhop, the potential volumic mass (Kg/m3) 
    242323      !! 
    243       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 
     324      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     325      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    244326      !!                Brown and Campana, Mon. Weather Rev., 1978 
    245327      !!---------------------------------------------------------------------- 
    … …  
    262344      SELECT CASE ( nn_eos ) 
    263345      ! 
    264       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
     346      CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==! 
    265347!CDIR NOVERRCHK 
    266348         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    … …  
    290372                  ! 
    291373                  ! add the compression terms 
     374                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     375                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     376                  zb = zbw + ze * zs 
     377                  ! 
     378                  zd = -1.480266e-4_wp 
     379                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
     380                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
     381                  za = ( zd*zsr + zc ) *zs + zaw 
     382                  ! 
     383                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     384                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     385                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
     386                  zk0= ( zb1*zsr + za1 )*zs + zkw 
     387                  ! 
     388                  ! masked in situ density anomaly 
     389                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
     390                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
     391               END DO 
     392            END DO 
     393         END DO 
     394         ! 
     395      CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==! 
     396!CDIR NOVERRCHK 
     397         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     398         ! 
     399         DO jk = 1, jpkm1 
     400            DO jj = 1, jpj 
     401               DO ji = 1, jpi 
     402                  zt = pts   (ji,jj,jk,jp_tem) 
     403                  zs = pts   (ji,jj,jk,jp_sal) 
     404                  zh = fsdept(ji,jj,jk)        ! depth 
     405                  zsr= zws   (ji,jj,jk)        ! square root salinity 
     406                  ! 
     407                  ! compute volumic mass pure water at atm pressure 
     408                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
     409                     &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
     410                  ! seawater volumic mass atm pressure 
     411                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   & 
     412                     &                                         -4.0899e-3_wp ) *zt+0.824493_wp 
     413                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     414                  zr4= 4.8314e-4_wp 
     415                  ! 
     416                  ! potential volumic mass (reference to the surface) 
     417                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     418                  ! 
     419                  ! save potential volumic mass 
     420                  prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk) 
     421                  ! 
     422                  ! add the compression terms 
    292423                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
    293424                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
    … …  
    323454         END DO 
    324455         ! 
     456      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     457         DO jk = 1, jpkm1 
     458            DO jj = 1, jpj 
     459               DO ji = 1, jpi 
     460                  zt = pts   (ji,jj,jk,jp_tem) - 10._wp 
     461                  zs = pts   (ji,jj,jk,jp_sal) - 35._wp 
     462                  zh = fsdept(ji,jj,jk)        ! depth 
     463                  ! masked in situ density anomaly 
     464                  prd(ji,jj,jk) = vallis_diff + vallis_ratio * (                                             & 
     465                     &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  & 
     466                     &   ) * tmask(ji,jj,jk) 
     467                  ! masked in situ density anomaly 
     468                 prhop(ji,jj,jk) = ( 1.0_wp - ( 1.67e-4_wp + 0.5e-5_wp*zt ) *zt + 0.78e-3_wp *zs )   & 
     469                     &              * 1027._wp * tmask(ji,jj,jk) 
     470                 ! 
     471               END DO 
     472            END DO 
     473         END DO 
     474         ! 
    325475      END SELECT 
    326476      ! 
    … …  
    342492      !!      defined through the namelist parameter nn_eos. * 2D field case 
    343493      !! 
    344       !! ** Method : 
    345       !!      nn_eos = 0 : Jackett and McDougall (1994) equation of state. 
     494      !! ** Method  :   5 cases: 
     495      !!      nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
     496      !!      nn_eos = 0 : standard NEMO equation of state, based on  
     497      !!         Jackett and McDougall (1995) equation of state. 
    346498      !!         the in situ density is computed directly as a function of 
    347499      !!         potential temperature relative to the surface (the opa t 
    … …  
    363515      !!               salinity 
    364516      !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 
     517      !!      nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 
     518      !!              prd(t,s,p) =  ( rho(t,s,p) - rau0 ) / rau0 
     519      !!              where the pressure p in decibars is approximated by the depth in meters. 
    365520      !!      Note that no boundary condition problem occurs in this routine 
    366521      !!      as pts are defined over the whole domain. 
    … …  
    368523      !! ** Action  : - prd , the in situ density (no units) 
    369524      !! 
    370       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 
     525      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     526      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
    371527      !!---------------------------------------------------------------------- 
    372528      !! 
    … …  
    391547      SELECT CASE( nn_eos ) 
    392548      ! 
    393       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
     549      CASE( -1 )                !==  Jackett and McDougall (1995) formulation  ==! 
    394550      ! 
    395551!CDIR NOVERRCHK 
    … …  
    403559            DO ji = 1, fs_jpim1   ! vector opt. 
    404560               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask 
    405                zt    = pts  (ji,jj,jp_tem)            ! interpolated T 
    406                zs    = pts  (ji,jj,jp_sal)            ! interpolated S 
     561               zt    = pts  (ji,jj,jp_tem)     ! interpolated T 
     562               zs    = pts  (ji,jj,jp_sal)     ! interpolated S 
     563               zsr   = zws  (ji,jj)            ! square root of interpolated S 
     564               zh    = pdep (ji,jj)            ! depth at the partial step level 
     565               ! 
     566               ! compute volumic mass pure water at atm pressure 
     567               zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   & 
     568                  &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 
     569               ! seawater volumic mass atm pressure 
     570               zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   & 
     571                  &                                   -4.0899e-3_wp ) *zt+0.824493_wp 
     572               zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 
     573               zr4 = 4.8314e-4_wp 
     574               ! 
     575               ! potential volumic mass (reference to the surface) 
     576               zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     577               ! 
     578               ! add the compression terms 
     579               ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     580               zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     581               zb = zbw + ze * zs 
     582               ! 
     583               zd = -1.480266e-4_wp 
     584               zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
     585               zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
     586               za = ( zd*zsr + zc ) *zs + zaw 
     587               ! 
     588               zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     589               za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     590               zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
     591               zk0= ( zb1*zsr + za1 )*zs + zkw 
     592               ! 
     593               ! masked in situ density anomaly 
     594               prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
     595            END DO 
     596         END DO 
     597         ! 
     598      CASE( 0 )                !==  modified Jackett and McDougall (1995) formulation  ==! 
     599      ! 
     600!CDIR NOVERRCHK 
     601         DO jj = 1, jpjm1 
     602!CDIR NOVERRCHK 
     603            DO ji = 1, fs_jpim1   ! vector opt. 
     604               zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 
     605            END DO 
     606         END DO 
     607         DO jj = 1, jpjm1 
     608            DO ji = 1, fs_jpim1   ! vector opt. 
     609               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask 
     610               zt    = pts  (ji,jj,jp_tem)     ! interpolated T 
     611               zs    = pts  (ji,jj,jp_sal)     ! interpolated S 
    407612               zsr   = zws  (ji,jj)            ! square root of interpolated S 
    408613               zh    = pdep (ji,jj)            ! depth at the partial step level 
    … …  
    455660         END DO 
    456661         ! 
     662      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     663         DO jj = 1, jpjm1 
     664            DO ji = 1, fs_jpim1   ! vector opt. 
     665               zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask 
     666               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! interpolated T 
     667               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! interpolated S 
     668               zh    = pdep (ji,jj)            ! depth at the partial step level 
     669               ! masked in situ density anomaly 
     670               prd(ji,jj) = vallis_diff + vallis_ratio * (                                                & 
     671                  &   -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh  & 
     672                  &   ) * zmask 
     673               ! 
     674            END DO 
     675         END DO 
     676         ! 
    457677      END SELECT 
    458678 
    … …  
    473693      !!      step of the input arguments 
    474694      !! 
    475       !! ** Method : 
    476       !!       * nn_eos = 0  : UNESCO sea water properties 
    477       !!         The brunt-vaisala frequency is computed using the polynomial 
    478       !!      polynomial expression of McDougall (1987): 
     695      !! ** Method  :   5 cases: 
     696      !!       * nn_eos = -1  : The brunt-vaisala frequency is computed for 
     697      !!       the Jackett and McDougall (1995) equation of state: 
    479698      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 
     699      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
     700      !!      computed and used in zdfddm module : 
     701      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 
     702      !!       * nn_eos = 0  : The brunt-vaisala frequency is based on the 
     703      !!            formulation of McDougall (1987). 
    480704      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
    481705      !!      computed and used in zdfddm module : 
    … …  
    494718      !! ** Action  : - pn2 : the brunt-vaisala frequency 
    495719      !! 
    496       !! References :   McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987. 
     720      !! References :   Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     721      !!                Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 
     722      !!                McDougall, J. Phys. Oceano., 1987 
    497723      !!---------------------------------------------------------------------- 
    498724      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius] 
    … …  
    501727      !! 
    502728      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    503       REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars 
     729      REAL(wp) ::   zgde3w, zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop      ! local scalars 
     730      REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM    !   -      - 
     731      REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom   !   -      - 
     732      REAL(wp) ::   zalpbet, zalpha, zbeta                                  !   -      - 
    504733#if defined key_zdfddm 
    505734      REAL(wp) ::   zds   ! local scalars 
    506735#endif 
     736      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    507737      !!---------------------------------------------------------------------- 
    508738 
    509739      ! 
    510740      IF( nn_timing == 1 ) CALL timing_start('bn2') 
     741      ! 
     742      CALL wrk_alloc( jpi, jpj, jpk, zws ) 
    511743      ! 
    512744      ! pn2 : interior points only (2=< jk =< jpkm1 ) 
    … …  
    515747      SELECT CASE( nn_eos ) 
    516748      ! 
    517       CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
     749      CASE( -1 )                !==  Jackett and McDougall (1995)  ==! 
     750         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     751         DO jk = 2, jpkm1 
     752            DO jj = 1, jpj 
     753               DO ji = 1, jpi 
     754                  zgde3w = grav / fse3w(ji,jj,jk) 
     755                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )         ! potential temperature at w-pt 
     756                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) )         ! salinity at w-pt 
     757                  zsr = 0.5 * ( zws(ji,jj,jk) + zws(ji,jj,jk-1) )        ! square root of S at w-pt 
     758                  zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point 
     759                  ! 
     760                  ! compute volumic mass pure water at atm pressure 
     761                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     762                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
     763                  ! seawater volumic mass atm pressure 
     764                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     765                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     766                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     767                  zr4= 4.8314e-4_wp 
     768                  ! 
     769                  ! potential volumic mass (reference to the surface) 
     770                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     771                  ! 
     772                  ! add the compression terms 
     773                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     774                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     775                  zb = zbw + ze * zs 
     776                  ! 
     777                  zd = -1.480266e-4_wp 
     778                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
     779                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
     780                  za = ( zd*zsr + zc ) *zs + zaw 
     781                  ! 
     782                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     783                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     784                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
     785                  zk0= ( zb1*zsr + za1 )*zs + zkw 
     786                  ! 
     787                  zM= ( zb*zh - za )*zh + zk0 
     788                  zDenom= 1._wp - zh / zM 
     789                  ! compute in-situ density 
     790                  zrho = zrhop/zDenom 
     791              ! alpha 
     792                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     793                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
     794                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     795                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     796                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     797                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     798                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     799                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     800                  zdDenom = zh * zdM / (zM*zM) 
     801                  zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 
     802                  ! beta 
     803                  zdzb  = ze 
     804                  zdza  = -3.0644505e-2_wp*zsr+zc 
     805                  zdzk0 = 1.5_wp*zb1*zsr+za1 
     806                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
     807                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     808                  zdDenom = zh * zdM / (zM*zM) 
     809                  zbeta =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 
     810                  ! alpha/beta 
     811                  zalpbet = zalpha / zbeta 
     812                  ! N2 
     813                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
     814                     &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
     815                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
     816#if defined key_zdfddm 
     817                  !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!! 
     818                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
     819                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
     820                  rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     821#endif 
     822               END DO 
     823            END DO 
     824         END DO 
     825         ! 
     826      CASE( 0 )                !==  McDougall (1987) formulation  ==! 
    518827         DO jk = 2, jpkm1 
    519828            DO jj = 1, jpj 
    … …  
    524833                  zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point 
    525834                  ! 
    526                   zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &   ! ratio alpha/beta 
     835                  zalpbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt  &   ! ratio alpha/beta 
    527836                     &                                  - 0.203814e-03_wp ) * zt   & 
    528837                     &                                  + 0.170907e-01_wp ) * zt   & 
    … …  
    550859                     ! 
    551860                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
    552                      &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
     861                     &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    553862                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
    554863#if defined key_zdfddm 
    … …  
    556865                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    557866                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
    558                   rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     867                  rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    559868#endif 
    560869               END DO 
    … …  
    564873      CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    565874         DO jk = 2, jpkm1 
    566             pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
     875            pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )        & 
     876              &                / fse3w(:,:,jk) * tmask(:,:,jk) 
    567877         END DO 
    568878         ! 
    … …  
    579889                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 
    580890                  IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 
    581                   rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     891                  rrau(ji,jj,jk) = rn_alpha / rn_beta * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
    582892               END DO 
    583893            END DO 
    584894         END DO 
    585895#endif 
     896         ! 
     897      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     898         DO jk = 2, jpkm1 
     899            DO jj = 1, jpj 
     900               DO ji = 1, jpi 
     901                  zgde3w = grav / fse3w(ji,jj,jk) 
     902                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) - 10._wp ! potential temperature at w-pt 
     903                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35._wp ! salinity anomaly (s-35) at w-pt 
     904                  zh = fsdepw(ji,jj,jk)                                                 ! depth in meters  at w-point 
     905                  ! 
     906                  zalpha = 1027._wp * ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) / rau0 ! alpha 
     907                  zbeta  = 1027._wp * 0.78e-3_wp / rau0  ! beta 
     908                  zalpbet = zalpha / zbeta  ! ratio alpha/beta 
     909                  ! 
     910                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
     911                     &          * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
     912                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
     913#if defined key_zdfddm 
     914                  !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!! 
     915                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
     916                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 
     917                  rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     918#endif 
     919               END DO 
     920            END DO 
     921         END DO 
     922         ! 
    586923      END SELECT 
    587924 
    … …  
    591928#endif 
    592929      ! 
     930      CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
     931      ! 
    593932      IF( nn_timing == 1 ) CALL timing_stop('bn2') 
    594933      ! 
    … …  
    603942      !! 
    604943      !! ** Method  :   calculates alpha / beta ratio at T-points 
    605       !!       * nn_eos = 0  : UNESCO sea water properties 
    606       !!                       The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 
    607       !!                       polynomial expression of McDougall (1987). 
     944      !!       * nn_eos = -1  : alpha / beta ratio is computed for 
     945      !!       the Jackett and McDougall (1995) equation of state: 
     946      !!                       Scalar beta0 is returned = 1. 
     947      !!       * nn_eos = 0  : alpha / beta ratio using the formulation of McDougall (1987). 
    608948      !!                       Scalar beta0 is returned = 1. 
    609949      !!       * nn_eos = 1  : linear equation of state (temperature only) 
    … …  
    611951      !!                       Scalar beta0 is returned = 0. 
    612952      !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    613       !!                       The alpha/beta ratio is returned as ralpbet 
     953      !!                       The alpha/beta ratio is returned as palpbet 
     954      !!                       Scalar beta0 is returned = 1. 
     955      !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006) 
    614956      !!                       Scalar beta0 is returned = 1. 
    615957      !! 
    … …  
    622964      !! 
    623965      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    624       REAL(wp) ::   zt, zs, zh   ! local scalars 
     966      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars 
     967      REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         - 
     968      REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 
     969      REAL(wp) ::   zalpha, zbeta                                ! local scalars 
     970      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
    625971      !!---------------------------------------------------------------------- 
    626972      ! 
    627973      IF( nn_timing == 1 ) CALL timing_start('eos_alpbet') 
    628974      ! 
     975      CALL wrk_alloc( jpi, jpj, jpk, zws ) 
     976      ! 
    629977      SELECT CASE ( nn_eos ) 
    630978      ! 
    631       CASE ( 0 )               ! Jackett and McDougall (1994) formulation 
     979      CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
     980         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     981         DO jk = 1, jpkm1 
     982            DO jj = 1, jpj 
     983               DO ji = 1, jpi 
     984                  zt = pts   (ji,jj,jk,jp_tem) 
     985                  zs = pts   (ji,jj,jk,jp_sal) 
     986                  zh = fsdept(ji,jj,jk)        ! depth 
     987                  zsr= zws   (ji,jj,jk)        ! square root salinity 
     988                  ! 
     989                  ! compute volumic mass pure water at atm pressure 
     990                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     991                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
     992                  ! seawater volumic mass atm pressure 
     993                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     994                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     995                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     996                  zr4= 4.8314e-4_wp 
     997                  ! 
     998                  ! potential volumic mass (reference to the surface) 
     999                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     1000                  ! 
     1001                  ! add the compression terms 
     1002                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     1003                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     1004                  zb = zbw + ze * zs 
     1005                  ! 
     1006                  zd = -1.480266e-4_wp 
     1007                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
     1008                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
     1009                  za = ( zd*zsr + zc ) *zs + zaw 
     1010                  ! 
     1011                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     1012                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     1013                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
     1014                  zk0= ( zb1*zsr + za1 )*zs + zkw 
     1015                  ! 
     1016                  zM= ( zb*zh - za )*zh + zk0 
     1017                  zDenom= 1._wp - zh / zM 
     1018                  ! compute in-situ density 
     1019                  zrho = zrhop/zDenom 
     1020              ! alpha 
     1021                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     1022                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
     1023                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     1024                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     1025                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     1026                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     1027                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     1028                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1029                  zdDenom = zh * zdM / (zM*zM) 
     1030                  zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 
     1031                  ! beta 
     1032                  zdzb  = ze 
     1033                  zdza  = -3.0644505e-2_wp*zsr+zc 
     1034                  zdzk0 = 1.5_wp*zb1*zsr+za1 
     1035                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
     1036                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1037                  zdDenom = zh * zdM / (zM*zM) 
     1038                  zbeta =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 
     1039                  ! alpha/beta 
     1040                  palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk) 
     1041               END DO 
     1042            END DO 
     1043         END DO 
     1044         beta0 = 1._wp 
     1045         ! 
     1046      CASE ( 0 )               ! McDougall (1987) formulation 
    6321047         DO jk = 1, jpk 
    6331048            DO jj = 1, jpj 
    … …  
    6601075         ! 
    6611076      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
    662          palpbet(:,:,:) = ralpbet 
     1077         palpbet(:,:,:) = rn_alpha / rn_beta 
     1078         beta0 = 1._wp 
     1079         ! 
     1080      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     1081         DO jk = 1, jpkm1 
     1082            DO jj = 1, jpj 
     1083               DO ji = 1, jpi 
     1084                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature 
     1085                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
     1086                  zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point 
     1087                  ! 
     1088                  zalpha = ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) ! alpha/vallis_ratio 
     1089                  zbeta  = 0.78e-3_wp  ! beta/vallis_ratio 
     1090                  ! 
     1091                  palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk) 
     1092               END DO 
     1093            END DO 
     1094         END DO 
    6631095         beta0 = 1._wp 
    6641096         ! 
    … …  
    6701102      END SELECT 
    6711103      ! 
     1104      CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
     1105      ! 
    6721106      IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet') 
    6731107      ! 
    … …  
    6751109 
    6761110 
     1111   SUBROUTINE eos_alpha_beta( pts, palpha, pbeta ) 
     1112      !!---------------------------------------------------------------------- 
     1113      !!                 ***  ROUTINE eos_alpha_beta  *** 
     1114      !! 
     1115      !! ** Purpose :   Calculates the in situ thermal/haline expansion terms at T-points 
     1116      !! 
     1117      !! ** Method  :   calculates alpha and beta at T-points 
     1118      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state 
     1119      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state 
     1120      !!       * nn_eos = 1  : linear equation of state (temperature only) 
     1121      !!                       palpha = rn_alpha 
     1122      !!                       pbeta = 0 
     1123      !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
     1124      !!                       palpha = rn_alpha 
     1125      !!                       pbeta = rn_beta 
     1126      !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006) 
     1127      !!                       alpha and beta ratios are returned as 3-D arrays. 
     1128      !! 
     1129      !! ** Action  : - palpha : thermal expansion ratio at T-points 
     1130      !!            : - pbeta  : haline expansion ratio at T-points 
     1131      !!---------------------------------------------------------------------- 
     1132      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts      ! pot. temperature & salinity 
     1133      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpha   ! alpha 
     1134      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pbeta    ! beta 
     1135      ! 
     1136      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     1137      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars 
     1138      REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         - 
     1139      REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 
     1140      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
     1141      !!---------------------------------------------------------------------- 
     1142      ! 
     1143      IF ( nn_timing == 1 )   CALL timing_start('eos_alpha_beta') 
     1144      ! 
     1145      CALL wrk_alloc( jpi, jpj, jpk, zws ) 
     1146      ! 
     1147      SELECT CASE ( nn_eos ) 
     1148      ! 
     1149      CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
     1150         ! 
     1151         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     1152         DO jk = 1, jpkm1 
     1153            DO jj = 1, jpj 
     1154               DO ji = 1, jpi 
     1155                  zt = pts   (ji,jj,jk,jp_tem) 
     1156                  zs = pts   (ji,jj,jk,jp_sal) 
     1157                  zh = fsdept(ji,jj,jk)        ! depth 
     1158                  zsr= zws   (ji,jj,jk)        ! square root salinity 
     1159                  ! 
     1160                  ! compute volumic mass pure water at atm pressure 
     1161                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     1162                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
     1163                  ! seawater volumic mass atm pressure 
     1164                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     1165                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     1166                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     1167                  zr4= 4.8314e-4_wp 
     1168                  ! 
     1169                  ! potential volumic mass (reference to the surface) 
     1170                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     1171                  ! 
     1172                  ! add the compression terms 
     1173                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     1174                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     1175                  zb = zbw + ze * zs 
     1176                  ! 
     1177                  zd = -1.480266e-4_wp 
     1178                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
     1179                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
     1180                  za = ( zd*zsr + zc ) *zs + zaw 
     1181                  ! 
     1182                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     1183                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     1184                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
     1185                  zk0= ( zb1*zsr + za1 )*zs + zkw 
     1186                  ! 
     1187                  zM= ( zb*zh - za )*zh + zk0 
     1188                  zDenom= 1._wp - zh / zM 
     1189                  ! compute in-situ density 
     1190                  zrho = zrhop/zDenom 
     1191              ! alpha 
     1192                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     1193                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
     1194                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     1195                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     1196                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     1197                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     1198                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     1199                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1200                  zdDenom = zh * zdM / (zM*zM) 
     1201                  palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 
     1202                  ! beta 
     1203                  zdzb  = ze 
     1204                  zdza  = -3.0644505e-2_wp*zsr+zc 
     1205                  zdzk0 = 1.5_wp*zb1*zsr+za1 
     1206                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
     1207                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1208                  zdDenom = zh * zdM / (zM*zM) 
     1209                  pbeta(ji,jj,jk) =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 
     1210               END DO 
     1211            END DO 
     1212         END DO 
     1213         ! 
     1214      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
     1215         ! 
     1216         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     1217         DO jk = 1, jpkm1 
     1218            DO jj = 1, jpj 
     1219               DO ji = 1, jpi 
     1220                  zt = pts   (ji,jj,jk,jp_tem) 
     1221                  zs = pts   (ji,jj,jk,jp_sal) 
     1222                  zh = fsdept(ji,jj,jk)        ! depth 
     1223                  zsr= zws   (ji,jj,jk)        ! square root salinity 
     1224                  ! 
     1225                  ! compute volumic mass pure water at atm pressure 
     1226                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     1227                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
     1228                  ! seawater volumic mass atm pressure 
     1229                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     1230                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     1231                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     1232                  zr4= 4.8314e-4_wp 
     1233                  ! 
     1234                  ! potential volumic mass (reference to the surface) 
     1235                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     1236                  ! 
     1237                  ! add the compression terms 
     1238                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     1239                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
     1240                  zb = zbw + ze * zs 
     1241                  ! 
     1242                  zd = -2.042967e-2_wp 
     1243                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
     1244                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
     1245                  za = ( zd*zsr + zc ) *zs + zaw 
     1246                  ! 
     1247                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     1248                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     1249                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
     1250                  zk0= ( zb1*zsr + za1 )*zs + zkw 
     1251                  ! 
     1252                  zM= ( zb*zh - za )*zh + zk0 
     1253                  zDenom= 1._wp - zh / zM 
     1254                  ! compute in-situ density 
     1255                  zrho = zrhop/zDenom 
     1256              ! alpha 
     1257                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     1258                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
     1259                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     1260                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     1261                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     1262                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     1263                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     1264                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1265                  zdDenom = zh * zdM / (zM*zM) 
     1266                  palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 
     1267                  ! beta 
     1268                  zdzb  = ze 
     1269                  zdza  = -3.0644505e-2_wp*zsr+zc 
     1270                  zdzk0 = 1.5_wp*zb1*zsr+za1 
     1271                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
     1272                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1273                  zdDenom = zh * zdM / (zM*zM) 
     1274                  pbeta(ji,jj,jk) =   ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 
     1275               END DO 
     1276            END DO 
     1277         END DO 
     1278         ! 
     1279      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==! 
     1280         palpha(:,:,:) = rn_alpha 
     1281         pbeta (:,:,:) = 0._wp 
     1282         ! 
     1283      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
     1284         palpha(:,:,:) = rn_alpha 
     1285         pbeta (:,:,:) = rn_beta 
     1286         ! 
     1287      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     1288         DO jk = 1, jpkm1 
     1289            DO jj = 1, jpj 
     1290               DO ji = 1, jpi 
     1291                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-10) 
     1292                  zh = fsdept(ji,jj,jk)                                                ! depth in meters  at w-point 
     1293                  palpha(ji,jj,jk) = vallis_ratio * & 
     1294                      &  ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk)   ! alpha 
     1295                  ! 
     1296               END DO 
     1297            END DO 
     1298         END DO 
     1299         pbeta (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:)   ! beta 
     1300         ! 
     1301      CASE DEFAULT 
     1302         IF(lwp) WRITE(numout,cform_err) 
     1303         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     1304         nstop = nstop + 1 
     1305         ! 
     1306      END SELECT 
     1307      ! 
     1308      CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
     1309      ! 
     1310      IF( nn_timing == 1 ) CALL timing_stop('eos_alpha_beta') 
     1311      ! 
     1312   END SUBROUTINE eos_alpha_beta 
     1313 
     1314 
     1315 
     1316   SUBROUTINE eos_pen( pts, palphaPE, pbetaPE, pPEanom ) 
     1317      !!---------------------------------------------------------------------- 
     1318      !!                 ***  ROUTINE eos_pen  *** 
     1319      !! 
     1320      !! ** Purpose :   Calculates alpha_PE, beta_PE and PE at T-points 
     1321      !! 
     1322      !! ** Method  :   PE is defined analytically as the vertical  
     1323      !!                   primitive of EOS times -g integrated between 0 and z>0. 
     1324      !!                PEanom is the PE anomaly: PEanom = ( PE - rau0 gz ) / rau0 gz 
     1325      !!                                           = -1/z * /int_z^0 rd , where rd density anomaly 
     1326      !!                dalphaPE and dbetaPE are partial derivatives of PE anomaly with respect to T and S: 
     1327      !!                    alphaPE = - 1/(rau0 gz) * dPE/dT = - dPEanom/dT 
     1328      !!                    betaPE  =   1/(rau0 gz) * dPE/dS = - dPEanom/dS 
     1329      !! 
     1330      !!       * nn_eos = -1 : Jackett and McDougall (1995) equation of state. 
     1331      !!       * nn_eos = 0  : modified Jackett and McDougall (1995) equation of state. 
     1332      !!       * nn_eos = 1  : linear equation of state (temperature only) 
     1333      !!                       palpha = rau0*g*rn_alpha*z 
     1334      !!                       pbeta = 0 
     1335      !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
     1336      !!                       palpha = rau0*g*rn_alpha*z 
     1337      !!                       pbeta = -rau0*g*rn_beta*z 
     1338      !!       * nn_eos = 3  : Vallis equation of state (Vallis 2006) 
     1339      !! 
     1340      !! ** Action  : - pPEanom     : PE anomaly given at T-points 
     1341      !!            : - palphaPE    : given at T-points 
     1342      !!            : - pbetaPE     : given at T-points 
     1343      !!---------------------------------------------------------------------- 
     1344      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity 
     1345      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palphaPE  ! partial derivative of PE anom 
     1346      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pbetaPE   ! with respect to T and S, resp. 
     1347      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pPEanom   ! potential energy anomaly 
     1348      ! 
     1349      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     1350      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop            ! temporary scalars 
     1351      REAL(wp) ::   ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM  !    -         - 
     1352      REAL(wp) ::   zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 
     1353      REAL(wp) ::   zap1, zdelta, zsdelta, zAp, zBp, zlogBp, zCp, zatanCp ! 
     1354      REAL(wp) ::   zddelta, zdAp, zdBp, zdlogBp, zdCp, zP, zdP, zEp      ! 
     1355      REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 
     1356      !!---------------------------------------------------------------------- 
     1357      ! 
     1358      IF ( nn_timing == 1 )   CALL timing_start('eos_pen') 
     1359      ! 
     1360      CALL wrk_alloc( jpi, jpj, jpk, zws ) 
     1361      ! 
     1362      SELECT CASE ( nn_eos ) 
     1363      ! 
     1364      CASE ( -1 )               ! Jackett and McDougall (1995) formulation 
     1365         ! 
     1366         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     1367         DO jk = 1, jpkm1 
     1368            DO jj = 1, jpj 
     1369               DO ji = 1, jpi 
     1370                  zt = pts   (ji,jj,jk,jp_tem) 
     1371                  zs = pts   (ji,jj,jk,jp_sal) 
     1372                  zh = fsdept(ji,jj,jk)        ! depth 
     1373                  zsr= zws   (ji,jj,jk)        ! square root salinity 
     1374                  ! 
     1375                  ! compute volumic mass pure water at atm pressure 
     1376                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     1377                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
     1378                  ! seawater volumic mass atm pressure 
     1379                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     1380                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     1381                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     1382                  zr4= 4.8314e-4_wp 
     1383                  ! 
     1384                  ! potential volumic mass (reference to the surface) 
     1385                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     1386                  ! 
     1387                  ! add the compression terms 
     1388                  ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 
     1389                  zbw= (  1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 
     1390                  zb = zbw + ze * zs 
     1391                  ! 
     1392                  zd = -1.480266e-4_wp 
     1393                  zc =   (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 
     1394                  zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 
     1395                  za = ( zd*zsr + zc ) *zs + zaw 
     1396                  ! 
     1397                  zb1=   (-4.619924e-3_wp*zt+9.085835e-2_wp )        *zt+3.886640_wp 
     1398                  za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp)       *zt-3.101089_wp ) *zt+528.4855_wp 
     1399                  zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 
     1400                  zk0= ( zb1*zsr + za1 )*zs + zkw 
     1401                  ! 
     1402                  zM= ( zb*zh - za )*zh + zk0 
     1403                  zDenom= 1._wp - zh / zM 
     1404                  ! compute in-situ density 
     1405                  zrho = zrhop/zDenom 
     1406                  ! 
     1407                  zap1    = 1._wp + za 
     1408                  zdelta  = 4._wp*zb*zk0 - zap1*zap1 
     1409                  zsdelta = sqrt( abs( zdelta ) ) 
     1410                  zAp     = zap1 / zb / zsdelta 
     1411                  zBp     = ( zM - zh ) / zk0 
     1412                  zlogBp  = log(abs(zBp)) 
     1413                  zCp     = zh * zsdelta / (2._wp*zk0-zh*zap1) 
     1414                  zatanCp = atan(zCp) 
     1415                  zP      = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb 
     1416                  ! compute potential energy 
     1417                  pPEanom(ji,jj,jk)     = ( ( zrhop * zP ) / rau0 / zh - 1._wp ) * tmask(ji,jj,jk) !!wrong 
     1418                  ! 
     1419              ! dPEdt 
     1420                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     1421                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
     1422                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     1423                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     1424                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     1425                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     1426                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     1427                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1428                  zdDenom = zh * zdM / (zM*zM) 
     1429                  ! 
     1430                  zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 
     1431                  zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 
     1432                  zdlogBp = zdM/(zM-zh) - zdzk0/zk0 
     1433                  zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 
     1434                  zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp)  
     1435                  ! 
     1436                  palphaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk)  
     1437                  ! 
     1438              ! dPEds 
     1439                  zdzb  = ze 
     1440                  zdza  = -3.0644505e-2_wp*zsr+zc 
     1441                  zdzk0 = 1.5_wp*zb1*zsr+za1 
     1442                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
     1443                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1444                  zdDenom = zh * zdM / (zM*zM) 
     1445                  ! 
     1446                  zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 
     1447                  zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 
     1448                  zdlogBp = zdM/(zM-zh) - zdzk0/zk0 
     1449                  zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 
     1450                  zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp)  
     1451                  ! 
     1452                  pbetaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 
     1453                  ! 
     1454               END DO 
     1455            END DO 
     1456         END DO 
     1457         ! 
     1458      CASE ( 0 )               ! modified Jackett and McDougall (1995) formulation 
     1459         ! 
     1460         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
     1461         DO jk = 1, jpkm1 
     1462            DO jj = 1, jpj 
     1463               DO ji = 1, jpi 
     1464                  zt = pts   (ji,jj,jk,jp_tem) 
     1465                  zs = pts   (ji,jj,jk,jp_sal) 
     1466                  zh = fsdept(ji,jj,jk)        ! depth 
     1467                  zsr= zws   (ji,jj,jk)        ! square root salinity 
     1468                  ! 
     1469                  ! compute volumic mass pure water at atm pressure 
     1470                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   & 
     1471                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp 
     1472                  ! seawater volumic mass atm pressure 
     1473                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        & 
     1474                     &                      -4.0899e-3_wp ) *zt+0.824493_wp 
     1475                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp 
     1476                  zr4= 4.8314e-4_wp 
     1477                  ! 
     1478                  ! potential volumic mass (reference to the surface) 
     1479                  zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 
     1480                  ! 
     1481                  ! add the compression terms 
     1482                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 
     1483                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 
     1484                  zb = zbw + ze * zs 
     1485                  ! 
     1486                  zd = -2.042967e-2_wp 
     1487                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 
     1488                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 
     1489                  za = ( zd*zsr + zc ) *zs + zaw 
     1490                  ! 
     1491                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp 
     1492                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp 
     1493                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 
     1494                  zk0= ( zb1*zsr + za1 )*zs + zkw 
     1495                  ! 
     1496                  zM= ( zb*zh - za )*zh + zk0 
     1497                  zDenom= 1._wp - zh / zM 
     1498                  ! compute in-situ density 
     1499                  zrho = zrhop/zDenom 
     1500                  ! 
     1501                  zap1    = 1._wp + za 
     1502                  zdelta  = 4._wp*zb*zk0 - zap1*zap1 
     1503                  zsdelta = sqrt( abs( zdelta ) ) 
     1504                  zAp     = zap1 / zb / zsdelta 
     1505                  zBp     = ( zM - zh ) / zk0 
     1506                  zlogBp  = log(abs(zBp)) 
     1507                  zCp     = zh * zsdelta / (2._wp*zk0-zh*zap1) 
     1508                  zatanCp = atan(zCp) 
     1509                  zP      = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb 
     1510                  ! compute potential energy 
     1511                  pPEanom(ji,jj,jk)     = - grav * zrhop * zP * tmask(ji,jj,jk) 
     1512                  ! 
     1513              ! dPEdt 
     1514                  zdzb  = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 
     1515                  zdza  = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 
     1516                  zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs  & 
     1517                     &            +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 
     1518                  zdrhop   = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr  & 
     1519                     &            +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs  & 
     1520                     &            +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 
     1521                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1522                  zdDenom = zh * zdM / (zM*zM) 
     1523                  ! 
     1524                  zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 
     1525                  zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 
     1526                  zdlogBp = zdM/(zM-zh) - zdzk0/zk0 
     1527                  zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 
     1528                  zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp)  
     1529                  ! 
     1530                  palphaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 
     1531                  ! 
     1532              ! dPEds 
     1533                  zdzb  = ze 
     1534                  zdza  = -3.0644505e-2_wp*zsr+zc 
     1535                  zdzk0 = 1.5_wp*zb1*zsr+za1 
     1536                  zdrhop   = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 
     1537                  zdM = (zdzb*zh-zdza)*zh+zdzk0 
     1538                  zdDenom = zh * zdM / (zM*zM) 
     1539                  ! 
     1540                  zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 
     1541                  zdAp    = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 
     1542                  zdlogBp = zdM/(zM-zh) - zdzk0/zk0 
     1543                  zdCp    = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 
     1544                  zdP     = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp)  
     1545                  ! 
     1546                  pbetaPE(ji,jj,jk)    = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 
     1547                  ! 
     1548               END DO 
     1549            END DO 
     1550         END DO 
     1551         ! 
     1552      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==! 
     1553         palphaPE(:,:,:) = rn_alpha * tmask(:,:,:) 
     1554         pbetaPE (:,:,:) = 0._wp 
     1555         DO jk = 1, jpkm1 
     1556            pPEanom(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
     1557         END DO 
     1558         ! 
     1559      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
     1560         palphaPE(:,:,:) = rn_alpha * tmask(:,:,:) 
     1561         pbetaPE (:,:,:) = rn_beta  * tmask(:,:,:) 
     1562         DO jk = 1, jpkm1 
     1563            pPEanom(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 
     1564         END DO 
     1565         ! 
     1566      CASE( 3 )                !==  Vallis (2006) simplified EOS  ==! 
     1567         DO jk = 1, jpkm1 
     1568            DO jj = 1, jpj 
     1569               DO ji = 1, jpi 
     1570                  zt = pts(ji,jj,jk,jp_tem) - 10._wp  ! potential temperature (t-10) 
     1571                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
     1572                  zh = fsdept(ji,jj,jk)                                                ! depth in meters  at w-point 
     1573                  ! 
     1574                  palphaPE(ji,jj,jk) = vallis_ratio * & 
     1575                      &  ( 1.67e-4_wp * ( 1._wp + 0.55e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk)   ! alphaPE 
     1576                  pPEanom(ji,jj,jk)     = vallis_diff + vallis_ratio *   & 
     1577                      &   ( - ( 1.67e-4_wp*(1._wp+0.55e-4_wp*zh) + 0.5e-5_wp*zt )*zt + 0.78e-3_wp*zs + 2.195e-6_wp*zh ) & 
     1578                      &   * tmask(ji,jj,jk) 
     1579                  ! 
     1580               END DO 
     1581            END DO 
     1582         END DO 
     1583         pbetaPE (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:)   ! betaPE=beta 
     1584         ! 
     1585      CASE DEFAULT 
     1586         IF(lwp) WRITE(numout,cform_err) 
     1587         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     1588         nstop = nstop + 1 
     1589         ! 
     1590      END SELECT 
     1591      ! 
     1592      CALL wrk_dealloc( jpi, jpj, jpk, zws ) 
     1593      ! 
     1594      IF( nn_timing == 1 ) CALL timing_stop('eos_pen') 
     1595      ! 
     1596   END SUBROUTINE eos_pen 
     1597 
     1598 
     1599 
    6771600   FUNCTION tfreez( psal ) RESULT( ptf ) 
    6781601      !!---------------------------------------------------------------------- 
    679       !!                 ***  ROUTINE eos_init  *** 
     1602      !!                 ***  ROUTINE tfreez  *** 
    6801603      !! 
    6811604      !! ** Purpose :   Compute the sea surface freezing temperature [Celcius] 
    … …  
    7241647      SELECT CASE( nn_eos )         ! check option 
    7251648      ! 
    726       CASE( 0 )                        !==  Jackett and McDougall (1994) formulation  ==! 
     1649      CASE( -1 )                        !==  Jackett and McDougall (1995) formulation  ==! 
    7271650         IF(lwp) WRITE(numout,*) 
    728          IF(lwp) WRITE(numout,*) '          use of Jackett & McDougall (1994) equation of state and' 
    729          IF(lwp) WRITE(numout,*) '                 McDougall (1987) Brunt-Vaisala frequency' 
     1651         IF(lwp) WRITE(numout,*) '          use of Jackett & McDougall (1995) equation of state' 
     1652         ! 
     1653      CASE( 0 )                        !==  modified Jackett and McDougall (1995) formulation  ==! 
     1654         IF(lwp) WRITE(numout,*) 
     1655         IF(lwp) WRITE(numout,*) '          use of modified Jackett & McDougall (1995) equation of state' 
     1656         IF(lwp) WRITE(numout,*) '                and McDougall (1987) Brunt-Vaisala frequency' 
    7301657         ! 
    7311658      CASE( 1 )                        !==  Linear formulation = F( temperature )  ==! 
    … …  
    7361663         ! 
    7371664      CASE( 2 )                        !==  Linear formulation = F( temperature , salinity )  ==! 
    738          ralpbet = rn_alpha / rn_beta 
    7391665         IF(lwp) WRITE(numout,*) 
    7401666         IF(lwp) WRITE(numout,*) '          use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )' 
     1667         ! 
     1668      CASE( 3 )                        !==  Vallis (2006) formulation  ==! 
     1669         IF(lwp) WRITE(numout,*) 
     1670         IF(lwp) WRITE(numout,*) '          use of Vallis (2006) equation of state' 
     1671         vallis_ratio = 1027._wp / rau0 
     1672         vallis_diff = (1027._wp-rau0) / rau0 
    7411673         ! 
    7421674      CASE DEFAULT                     !==  ERROR in nn_eos  ==! 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90

    r3718 r3876  
    44   !! Ocean  tracers:  horizontal & vertical advective trend 
    55   !!====================================================================== 
    6    !! History :  8.2  ! 2001-08  (G. Madec, E. Durand) trahad+trazad=traadv  
    7    !!            1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
    8    !!            9.0  ! 2004-08  (C. Talandier) New trends organization 
     6   !! History :  OPA  ! 2001-08  (G. Madec, E. Durand) v8.2 trahad+trazad=traadv  
     7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
     8   !!             -   ! 2004-08  (C. Talandier) New trends organization 
    99   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization 
    1010   !!            2.0  ! 2006-04  (R. Benshila, G. Madec) Step reorganization 
    … …  
    2121   USE dom_oce         ! ocean space and time domain 
    2222   USE eosbn2          ! equation of state 
    23    USE trdmod_oce      ! tracers trends 
    24    USE trdtra          ! tracers trends 
     23   USE trd_oce         ! trends: ocean variables 
     24   USE trdtra          ! trends manager: tracers  
    2525   USE closea          ! closed sea 
    2626   USE sbcrnf          ! river runoffs 
    … …  
    3737   PRIVATE 
    3838 
    39    PUBLIC   tra_adv_cen2       ! routine called by step.F90 
    40    PUBLIC   ups_orca_set       ! routine used by traadv_cen2_jki.F90 
     39   PUBLIC   tra_adv_cen2   ! routine called by step.F90 
     40   PUBLIC   ups_orca_set   ! routine used by traadv_cen2_jki.F90 
    4141 
    4242   LOGICAL  :: l_trd       ! flag to compute trends 
    … …  
    5555 
    5656   SUBROUTINE tra_adv_cen2( kt, kit000, cdtype, pun, pvn, pwn,     & 
    57       &                                 ptb, ptn, pta, kjpt   )  
     57      &                                         ptb, ptn, pta, kjpt   )  
    5858      !!---------------------------------------------------------------------- 
    5959      !!                  ***  ROUTINE tra_adv_cen2  *** 
    … …  
    8585      !!       * Add this trend now to the general trend of tracer (ta,sa): 
    8686      !!               pta = pta + ztra 
    87       !!       * trend diagnostic ('key_trdtra' defined): the trend is 
     87      !!       * trend diagnostic (l_trdtra=T or l_trctra=T): the trend is 
    8888      !!      saved for diagnostics. The trends saved is expressed as 
    89       !!      Uh.gradh(T), i.e. 
    90       !!                     save trend = ztra + ptn divn 
     89      !!      Uh.gradh(T), i.e.  save trend = ztra + ptn divn 
    9190      !! 
    9291      !!         Part II : vertical advection 
    … …  
    104103      !!         Add this trend now to the general trend of tracer (ta,sa): 
    105104      !!             pta = pta + ztra 
    106       !!         Trend diagnostic ('key_trdtra' defined): the trend is 
     105      !!         Trend diagnostic (l_trdtra=T or l_trctra=T): the trend is 
    107106      !!      saved for diagnostics. The trends saved is expressed as : 
    108107      !!             save trend =  w.gradz(T) = ztra - ptn divn. 
    … …  
    144143         IF(lwp) WRITE(numout,*) 
    145144         ! 
    146          IF ( .NOT. ALLOCATED( upsmsk ) )  THEN 
     145         IF( .NOT. ALLOCATED( upsmsk ) )  THEN 
    147146             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 
    148147             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array') 
    … …  
    260259         END DO 
    261260 
    262          !                                 ! trend diagnostics (contribution of upstream fluxes) 
     261         !                                 ! trend diagnostics 
    263262         IF( l_trd ) THEN 
    264             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 
    265             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    266             CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) ) 
     263            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
     264            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     265            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
    267266         END IF 
    268267         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    … …  
    272271         ENDIF 
    273272         ! 
    274       ENDDO 
     273      END DO 
    275274 
    276275      ! ---------------------------  required in restart file to ensure restartability) 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

    r3718 r3876  
    1616   !!---------------------------------------------------------------------- 
    1717   USE oce            ! ocean dynamics and active tracers 
     18   USE trc_oce        ! share passive tracers/Ocean variables 
    1819   USE dom_oce        ! ocean space and time domain 
    19    USE trdmod_oce     ! tracers trends  
    20    USE trdtra         ! tracers trends  
    21    USE in_out_manager ! I/O manager 
     20   USE trd_oce        ! trends: ocean variables 
     21   USE trdtra         ! tracers trends manager 
    2222   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    23    USE trabbl         ! tracers: bottom boundary layer 
    24    USE lib_mpp        ! distribued memory computing 
    25    USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
     23   USE sbcrnf          ! river runoffs 
    2624   USE diaptr         ! poleward transport diagnostics 
    27    USE trc_oce        ! share passive tracers/Ocean variables 
     25   ! 
    2826   USE wrk_nemo       ! Memory Allocation 
    2927   USE timing         ! Timing 
    3028   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    31    USE eosbn2          ! equation of state 
    32    USE sbcrnf          ! river runoffs 
     29   USE in_out_manager ! I/O manager 
     30   USE lib_mpp        ! distribued memory computing 
     31   USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
     32!!gm   USE trabbl         ! tracers: bottom boundary layer 
     33!!gm   USE eosbn2          ! equation of state 
    3334 
    3435   IMPLICIT NONE 
    … …  
    191192                  zalpha = 0.5 - z0u 
    192193                  zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 
    193                   zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * (zu * zslpx(ji+1,jj,jk)) 
    194                   zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * (zu * zslpx(ji  ,jj,jk)) 
     194                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 
     195                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk) 
    195196                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    196197                  ! 
    … …  
    198199                  zalpha = 0.5 - z0v 
    199200                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 
    200                   zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * (zv * zslpy(ji,jj+1,jk)) 
    201                   zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * (zv * zslpy(ji,jj  ,jk)) 
     201                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 
     202                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk) 
    202203                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    203204               END DO 
    … …  
    222223         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    223224         IF( l_trd )  THEN 
    224             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptb(:,:,:,jn) ) 
    225             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptb(:,:,:,jn) ) 
     225            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 
     226            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
    226227         END IF 
    227228         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    … …  
    273274                  zalpha = 0.5 + z0w 
    274275                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr  
    275                   zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * (zw * zslpx(ji,jj,jk+1)) 
    276                   zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * (zw * zslpx(ji,jj,jk  )) 
     276                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
     277                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
    277278                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    278279               END DO  
    … …  
    293294         END DO 
    294295         !                                 ! Save the vertical advective trends for diagnostic 
    295          IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwx, pwn, ptb(:,:,:,jn) ) 
     296         IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
    296297         ! 
    297298      ENDDO 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90

    r3625 r3876  
    1313   !!---------------------------------------------------------------------- 
    1414   USE oce             ! ocean dynamics and active tracers 
     15   USE trc_oce         ! share passive tracers/Ocean variables 
    1516   USE dom_oce         ! ocean space and time domain 
    16    USE trdmod_oce      ! tracers trends  
    17    USE trdtra          ! tracers trends  
     17   USE trd_oce         ! trends: ocean variables 
     18   USE trdtra          ! trends manager: tracers  
    1819   USE in_out_manager  ! I/O manager 
    1920   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient 
    20    USE trabbl          ! tracers: bottom boundary layer 
     21!!gm   USE trabbl          ! tracers: bottom boundary layer 
     22   USE diaptr          ! poleward transport diagnostics 
     23   ! 
    2124   USE lib_mpp         ! distribued memory computing 
    2225   USE lbclnk          ! ocean lateral boundary condition (or mpp link)  
    23    USE diaptr          ! poleward transport diagnostics 
    24    USE trc_oce         ! share passive tracers/Ocean variables 
    2526   USE wrk_nemo        ! Memory Allocation 
    2627   USE timing          ! Timing 
    … …  
    201202         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    202203         IF( l_trd ) THEN 
    203             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptb(:,:,:,jn) ) 
    204             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptb(:,:,:,jn) ) 
     204            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 
     205            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
    205206         END IF 
    206207 
    … …  
    284285         END DO 
    285286         !                       ! trend diagnostics (contribution of upstream fluxes) 
    286          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwx, pwn, ptb(:,:,:,jn) ) 
     287         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
    287288         ! 
    288289      END DO 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r3625 r3876  
    1717   USE oce             ! ocean dynamics and active tracers 
    1818   USE dom_oce         ! ocean space and time domain 
    19    USE trdmod_oce      ! ocean space and time domain 
    20    USE trdtra          ! ocean tracers trends  
    21    USE trabbl          ! advective term in the BBL 
     19   USE trc_oce         ! share passive tracers/Ocean variables 
     20   USE trd_oce         ! trends: ocean variables 
     21   USE trdtra          ! trends manager: tracers  
     22!!gm   USE trabbl          ! advective term in the BBL 
     23   USE dynspg_oce      ! surface pressure gradient variables 
     24   USE diaptr          ! poleward transport diagnostics 
     25   ! 
    2226   USE lib_mpp         ! distribued memory computing 
    2327   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    24    USE dynspg_oce      ! surface pressure gradient variables 
    2528   USE in_out_manager  ! I/O manager 
    26    USE diaptr          ! poleward transport diagnostics 
    27    USE trc_oce         ! share passive tracers/Ocean variables 
    2829   USE wrk_nemo        ! Memory Allocation 
    2930   USE timing          ! Timing 
    … …  
    233234         END DO 
    234235         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    235          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 
     236         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
    236237         ! 
    237238      END DO 
    … …  
    359360         END DO 
    360361         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    361          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     362         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    362363         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    363364         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    … …  
    422423         END DO 
    423424         !                                 ! Save the vertical advective trends for diagnostic 
    424          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) ) 
     425         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
    425426         ! 
    426427      END DO 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r3625 r3876  
    2222   USE oce            ! ocean dynamics and active tracers 
    2323   USE dom_oce        ! ocean space and time domain 
    24    USE trdmod_oce     ! tracers trends 
     24   USE trc_oce        ! share passive tracers/Ocean variables 
     25   USE trd_oce        ! trends: ocean variables 
    2526   USE trdtra         ! tracers trends 
    26    USE in_out_manager ! I/O manager 
    2727   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
     28   USE diaptr         ! poleward transport diagnostics 
     29   ! 
    2830   USE lib_mpp        ! MPP library 
    2931   USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    30    USE diaptr         ! poleward transport diagnostics 
    31    USE trc_oce        ! share passive tracers/Ocean variables 
     32   USE in_out_manager ! I/O manager 
    3233   USE wrk_nemo       ! Memory Allocation 
    3334   USE timing         ! Timing 
    … …  
    228229            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
    229230             
    230             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )    
    231             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )   
    232             CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )  
     231            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )    
     232            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )   
     233            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )  
    233234         END IF 
    234235         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r3787 r3876  
    1414   USE oce            ! ocean dynamics and active tracers 
    1515   USE dom_oce        ! ocean space and time domain 
    16    USE trdmod_oce     ! ocean space and time domain 
    17    USE trdtra 
    18    USE lib_mpp 
     16   USE trc_oce        ! share passive tracers/Ocean variables 
     17   USE trd_oce        ! trends: ocean variables 
     18   USE trdtra         ! trends manager: tracers  
     19   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
     20   USE diaptr         ! poleward transport diagnostics 
     21   ! 
     22   USE lib_mpp        ! I/O library 
    1923   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    2024   USE in_out_manager ! I/O manager 
    21    USE diaptr         ! poleward transport diagnostics 
    22    USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    23    USE trc_oce        ! share passive tracers/Ocean variables 
    2425   USE wrk_nemo       ! Memory Allocation 
    2526   USE timing         ! Timing 
    … …  
    5152      !!      and add it to the general trend of passive tracer equations. 
    5253      !! 
    53       !! ** Method  :   The upstream biased 3rd order scheme (UBS) is based on an 
     54      !! ** Method  :   The upstream biased scheme (UBS) is based on a 3rd order 
    5455      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 
    5556      !!      It is only used in the horizontal direction. 
    … …  
    182183         !                                 ! trend diagnostics (contribution of upstream fluxes) 
    183184         IF( l_trd ) THEN 
    184              CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 
    185              CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     185             CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
     186             CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    186187         END IF 
    187188         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    … …  
    265266               END DO 
    266267            END DO 
    267             CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv ) 
     268            CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv ) 
    268269         ENDIF 
    269270         ! 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

    r3625 r3876  
    1818   USE dom_oce         ! domain: ocean 
    1919   USE phycst          ! physical constants 
    20    USE trdmod_oce      ! trends: ocean variables  
    21    USE trdtra          ! trends: active tracers  
     20   USE trd_oce         ! trends: ocean variables 
     21   USE trdtra          ! trends manager: tracers  
    2222   USE in_out_manager  ! I/O manager 
    2323   USE prtctl          ! Print control 
    … …  
    9999      IF( l_trdtra ) THEN        ! Save the geothermal heat flux trend for diagnostics 
    100100         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    101          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbc, ztrdt ) 
     101         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 
    102102         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
    103103      ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r3764 r3876  
    2828   USE phycst         ! physical constant 
    2929   USE eosbn2         ! equation of state 
    30    USE trdmod_oce     ! trends: ocean variables 
    31    USE trdtra         ! trends: active tracers 
    32    USE iom            ! IOM server 
     30   USE trd_oce        ! trends: ocean variables 
     31   USE trdtra         ! trends manager: tracers  
     32   USE iom            ! IOM server                
    3333   USE in_out_manager ! I/O manager 
    3434   USE lbclnk         ! ocean lateral boundary conditions 
    … …  
    148148         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    149149         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    150          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt ) 
    151          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds ) 
    152          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
     150         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 
     151         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) 
     152         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )  
    153153      ENDIF 
    154154      ! 
    … …  
    180180      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 
    181181      !!---------------------------------------------------------------------- 
    182       ! 
    183182      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers 
    184183      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields 
    … …  
    399398                                          - 0.121555e-07 ) * pfh 
    400399      !!---------------------------------------------------------------------- 
    401  
    402400      ! 
    403401      IF( nn_timing == 1 )  CALL timing_start( 'bbl') 
    … …  
    405403      CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
    406404      ! 
    407  
    408405      IF( kt == kit000 )  THEN 
    409406         IF(lwp)  WRITE(numout,*) 
    … …  
    411408         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~' 
    412409      ENDIF 
    413  
    414410      !                                        !* bottom temperature, salinity, velocity and depth 
    415411#if defined key_vectopt_loop 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r3294 r3876  
    2727   USE oce            ! ocean: variables 
    2828   USE dom_oce        ! ocean: domain variables 
    29    USE trdmod_oce     ! ocean: trend variables 
    30    USE trdtra         ! active tracers: trends 
     29   USE trd_oce        ! trends: ocean variables 
     30   USE trdtra         ! trends manager: tracers  
    3131   USE zdf_oce        ! ocean: vertical physics 
    3232   USE phycst         ! physical constants 
    … …  
    4747   PUBLIC   dtacof_zoom  ! routine called by in both tradmp.F90 and trcdmp.F90 
    4848 
    49    !                                !!* Namelist namtra_dmp : T & S newtonian damping * 
     49   !                                         !!* Namelist namtra_dmp : T & S newtonian damping * 
    5050   LOGICAL, PUBLIC ::   ln_tradmp = .TRUE.    !: internal damping flag 
    5151   INTEGER         ::   nn_hdmp   =   -1      ! = 0/-1/'latitude' for damping over T and S 
    … …  
    111111      ! 
    112112      CALL wrk_alloc( jpi, jpj, jpk, jpts,  zts_dta ) 
     113      ! 
    113114      !                           !==   input T-S data at kt   ==! 
    114115      CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt 
    … …  
    171172      ! 
    172173      IF( l_trdtra )   THEN       ! trend diagnostic 
    173          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_dmp, ttrdmp ) 
    174          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_dmp, strdmp ) 
     174         CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ttrdmp ) 
     175         CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, strdmp ) 
    175176      ENDIF 
    176177      !                           ! Control print 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r3294 r3876  
    2323   USE traldf_iso_grif ! lateral mixing          (tra_ldf_iso_grif routine) 
    2424   USE traldf_lap      ! lateral mixing               (tra_ldf_lap routine) 
    25    USE trdmod_oce      ! ocean space and time domain 
    26    USE trdtra          ! ocean active tracers trends 
     25   USE trd_oce         ! trends: ocean variables 
     26   USE trdtra          ! trends manager: tracers  
    2727   USE prtctl          ! Print control 
    2828   USE in_out_manager  ! I/O manager 
    … …  
    3535   PRIVATE 
    3636 
    37    PUBLIC   tra_ldf         ! called by step.F90  
    38    PUBLIC   tra_ldf_init    ! called by opa.F90  
     37   PUBLIC   tra_ldf        ! called by step.F90  
     38   PUBLIC   tra_ldf_init   ! called by opa.F90  
    3939   ! 
    4040   INTEGER ::   nldf = 0   ! type of lateral diffusion used defined from ln_traldf_... namlist logicals) 
    … …  
    112112         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    113113         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    114          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_ldf, ztrdt ) 
    115          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_ldf, ztrds ) 
     114         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
     115         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 
    116116         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )  
    117117      ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

    r3294 r3876  
    1717   USE dom_oce         ! ocean space and time domain 
    1818   USE zdf_oce         ! ocean vertical physics 
    19    USE trdmod_oce      ! ocean active tracer trends 
    20    USE trdtra          ! ocean active tracer trends 
     19   USE trd_oce         ! trends: ocean variables 
     20   USE trdtra          ! trends manager: tracers  
    2121   USE eosbn2          ! equation of state (eos routine)  
    2222   USE lbclnk          ! lateral boundary conditions (or mpp link) 
    … …  
    5454      !! 
    5555      !! ** Action  : - (ta,sa) after the application od the npc scheme 
    56       !!              - save the associated trends (ttrd,strd) ('key_trdtra') 
     56      !!              - save the associated trends (l_trdtra=T) 
    5757      !! 
    5858      !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371. 
    … …  
    196196         !                                                ! =============== 
    197197         !  
     198       
     199         ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign) 
     200         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 
     201 
    198202         IF( l_trdtra )   THEN         ! save the Non penetrative mixing trends for diagnostic 
    199203            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    200204            ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    201             CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_npc, ztrdt ) 
    202             CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_npc, ztrds ) 
     205            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 
     206            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 
    203207            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    204208         ENDIF 
    205209       
    206          ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign) 
    207          ! ------------------------------============ 
    208          CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 
    209        
    210  
    211          !  2. non penetrative convective scheme statistics 
    212          !  ----------------------------------------------- 
     210         ! non penetrative convective scheme statistics 
    213211         IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN 
    214212            IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable',   & 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r3294 r3876  
    2727   USE dom_oce         ! ocean space and time domain variables  
    2828   USE sbc_oce         ! surface boundary condition: ocean 
    29    USE zdf_oce         ! ??? 
     29   USE zdf_oce         ! ocean vertical mixing 
    3030   USE domvvl          ! variable volume 
    3131   USE dynspg_oce      ! surface     pressure gradient variables 
    3232   USE dynhpg          ! hydrostatic pressure gradient  
    33    USE trdmod_oce      ! ocean space and time domain variables  
    34    USE trdtra          ! ocean active tracers trends  
    35    USE phycst 
    36    USE obc_oce 
     33   USE trd_oce         ! trends: ocean variables 
     34   USE trdtra          ! trends manager: tracers  
     35   USE traqsr          ! penetrative solar radiation (needed for nksr) 
     36   USE phycst          ! physical constant 
     37   USE ldftra_oce      ! lateral physics on tracers 
     38   USE bdy_oce         ! BDY open boundary condition variables 
     39   USE bdytra          ! BDY open boundary condition (bdy_tra routine) 
     40   USE obc_oce         ! open boundary condition variables 
    3741   USE obctra          ! open boundary condition (obc_tra routine) 
    38    USE bdy_oce 
    39    USE bdytra          ! open boundary condition (bdy_tra routine) 
     42   ! 
    4043   USE in_out_manager  ! I/O manager 
    4144   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    4245   USE prtctl          ! Print control 
    43    USE traqsr          ! penetrative solar radiation (needed for nksr) 
     46   USE wrk_nemo        ! Memory allocation 
     47   USE timing          ! Timing 
    4448#if defined key_agrif 
    4549   USE agrif_opa_update 
    4650   USE agrif_opa_interp 
    4751#endif 
    48    USE wrk_nemo        ! Memory allocation 
    49    USE timing          ! Timing 
    5052 
    5153   IMPLICIT NONE 
    … …  
    132134         ztrdt(:,:,:) = tsn(:,:,:,jp_tem)  
    133135         ztrds(:,:,:) = tsn(:,:,:,jp_sal) 
     136         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend  
     137            CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 
     138            CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds ) 
     139         ENDIF 
    134140      ENDIF 
    135141 
    … …  
    159165            ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact 
    160166         END DO 
    161          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt ) 
    162          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds ) 
     167         CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt ) 
     168         CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds ) 
    163169         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    164170      END IF 
    … …  
    168174         &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask ) 
    169175      ! 
    170       ! 
    171       IF( nn_timing == 1 )  CALL timing_stop('tra_nxt') 
     176      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt') 
    172177      ! 
    173178   END SUBROUTINE tra_nxt 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r3680 r3876  
    1313 
    1414   !!---------------------------------------------------------------------- 
    15    !!   tra_qsr      : trend due to the solar radiation penetration 
    16    !!   tra_qsr_init : solar radiation penetration initialization 
     15   !!   tra_qsr       : trend due to the solar radiation penetration 
     16   !!   tra_qsr_init  : solar radiation penetration initialization 
    1717   !!---------------------------------------------------------------------- 
    18    USE oce             ! ocean dynamics and active tracers 
    19    USE dom_oce         ! ocean space and time domain 
    20    USE sbc_oce         ! surface boundary condition: ocean 
    21    USE trc_oce         ! share SMS/Ocean variables 
    22    USE trdmod_oce      ! ocean variables trends 
    23    USE trdtra          ! ocean active tracers trends  
    24    USE in_out_manager  ! I/O manager 
    25    USE phycst          ! physical constants 
    26    USE prtctl          ! Print control 
    27    USE iom             ! I/O manager 
    28    USE fldread         ! read input fields 
    29    USE lib_mpp         ! MPP library 
     18   USE oce            ! ocean dynamics and active tracers 
     19   USE dom_oce        ! ocean space and time domain 
     20   USE sbc_oce        ! surface boundary condition: ocean 
     21   USE trc_oce        ! share SMS/Ocean variables 
     22   USE trd_oce        ! trends: ocean variables 
     23   USE trdtra         ! trends manager: tracers  
     24   USE phycst         ! physical constants 
     25   ! 
     26   USE in_out_manager ! I/O manager 
     27   USE prtctl         ! Print control 
     28   USE iom            ! I/O manager 
     29   USE fldread        ! read input fields 
     30   USE lib_mpp        ! MPP library 
    3031   USE wrk_nemo       ! Memory Allocation 
    3132   USE timing         ! Timing 
    32  
    3333 
    3434   IMPLICIT NONE 
    … …  
    4848   REAL(wp), PUBLIC ::   rn_si1     = 23.0_wp   !: deepest depth of extinction (water type I)       (2 bands) 
    4949    
    50    ! Module variables 
    51    REAL(wp) ::   xsi0r                           !: inverse of rn_si0 
    52    REAL(wp) ::   xsi1r                           !: inverse of rn_si1 
     50   INTEGER , PUBLIC ::   nksr   !: levels below which the light cannot penetrate ( depth larger than 391 m) 
     51 
     52   REAL(wp)                  ::   xsi0r, xsi1r        ! inverse of rn_si0 and rn_si1, resp. 
     53   REAL(wp), DIMENSION(3,61) ::   rkrgb               ! tabulated attenuation coefficients for RGB absorption 
    5354   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
    54    INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m) 
    55    REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption 
    5655 
    5756   !! * Substitutions 
    … …  
    8786      !! 
    8887      !! ** Action  : - update ta with the penetrative solar radiation trend 
    89       !!              - save the trend in ttrd ('key_trdtra') 
     88      !!              - send the trend to trdtra (l_trdtra=T) 
    9089      !! 
    9190      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 
    9291      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 
    9392      !!---------------------------------------------------------------------- 
    94       ! 
    9593      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
    9694      ! 
    … …  
    116114      ENDIF 
    117115 
    118       IF( l_trdtra ) THEN      ! Save ta and sa trends 
     116      IF( l_trdtra ) THEN      ! Save temperature trend 
    119117         CALL wrk_alloc( jpi, jpj, jpk, ztrdt )  
    120118         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    … …  
    141139      !                                        Compute now qsr tracer content field 
    142140      !                                        ************************************ 
    143        
    144141      !                                           ! ============================================== ! 
    145142      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  ! 
    … …  
    168165            IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume 
    169166               ! 
    170                IF( nn_chldta == 1 ) THEN                             !* Variable Chlorophyll 
     167               IF( nn_chldta == 1 ) THEN                             !- Variable Chlorophyll 
    171168                  ! 
    172169                  CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step 
    … …  
    184181                     END DO 
    185182                  END DO 
    186                ELSE                                            ! Variable ocean volume but constant chrlorophyll 
    187                   zchl = 0.05                                     ! constant chlorophyll 
     183               ELSE                                                  !- Variable ocean volume but constant chrlorophyll 
     184                  zchl = 0.05                                           ! constant chlorophyll 
    188185                  irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 
    189                   zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll  
     186                  zekb(:,:) = rkrgb(1,irgb)                             ! Separation in R-G-B depending of the chlorophyll  
    190187                  zekg(:,:) = rkrgb(2,irgb) 
    191188                  zekr(:,:) = rkrgb(3,irgb) 
    192189               ENDIF 
    193190               ! 
    194                zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B 
    195                ze0(:,:,1) = rn_abs  * qsr(:,:) 
    196                ze1(:,:,1) = zcoef * qsr(:,:) 
    197                ze2(:,:,1) = zcoef * qsr(:,:) 
    198                ze3(:,:,1) = zcoef * qsr(:,:) 
    199                zea(:,:,1) =         qsr(:,:) 
     191               zcoef  = ( 1. - rn_abs ) / 3.e0                       !- equi-partition in R-G-B 
     192               ze0(:,:,1) = rn_abs * qsr(:,:) 
     193               ze1(:,:,1) =  zcoef * qsr(:,:) 
     194               ze2(:,:,1) =  zcoef * qsr(:,:) 
     195               ze3(:,:,1) =  zcoef * qsr(:,:) 
     196               zea(:,:,1) =          qsr(:,:) 
    200197               ! 
    201198               DO jk = 2, nksr+1 
    … …  
    283280      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics 
    284281         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    285          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_qsr, ztrdt ) 
     282         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
    286283         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )  
    287284      ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r3764 r3876  
    1818   USE dom_oce         ! ocean space domain variables 
    1919   USE phycst          ! physical constant 
     20   USE sbcmod          ! ln_rnf   
     21   USE sbcrnf          ! River runoff   
    2022   USE traqsr          ! solar radiation penetration 
    21    USE trdmod_oce      ! ocean trends  
    22    USE trdtra          ! ocean trends 
     23   USE trd_oce         ! trends: ocean variables 
     24   USE trdtra          ! trends manager: tracers  
     25   ! 
    2326   USE in_out_manager  ! I/O manager 
    2427   USE prtctl          ! Print control 
    25    USE sbcrnf          ! River runoff   
    26    USE sbcmod          ! ln_rnf   
    27    USE iom 
     28   USE iom             ! I/O library 
    2829   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2930   USE wrk_nemo        ! Memory Allocation 
    … …  
    9192      !!         where emp, the surface freshwater budget (evaporation minus 
    9293      !!         precipitation minus runoff) given in kg/m2/s is divided 
    93       !!         by rau0 = 1020 kg/m3 (density of sea water) to obtain m/s.     
     94      !!         by rau0 = 1035 kg/m3 (density of sea water) to obtain m/s.     
    9495      !!         Note: even though Fwe does not appear explicitly for  
    9596      !!         temperature in this routine, the heat carried by the water 
    … …  
    107108      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated 
    108109      !!                with the tracer surface boundary condition  
    109       !!              - save the trend it in ttrd ('key_trdtra') 
     110      !!              - send trends to trdtra module (l_trdtra=T) 
    110111      !!---------------------------------------------------------------------- 
    111112      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    … …  
    124125      ENDIF 
    125126 
    126       IF( l_trdtra )   THEN                    !* Save ta and sa trends 
     127      IF( l_trdtra ) THEN                    !* Save ta and sa trends 
    127128         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )  
    128129         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    … …  
    137138 
    138139      !---------------------------------------- 
    139       !        EMP, EMPS and QNS effects 
     140      !        EMP, SFX and QNS effects 
    140141      !---------------------------------------- 
    141142      !                                          Set before sbc tracer content fields 
    … …  
    146147              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 
    147148            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file' 
    148             zfact = 0.5e0 
     149            zfact = 0.5_wp 
    149150            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend 
    150151            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend 
    151152         ELSE                                         ! No restart or restart not found: Euler forward time stepping 
    152             zfact = 1.e0 
    153             sbc_tsc_b(:,:,:) = 0.e0 
     153            zfact = 1._wp 
     154            sbc_tsc_b(:,:,:) = 0._wp 
    154155         ENDIF 
    155156      ELSE                                         ! Swap of forcing fields 
    156157         !                                         ! ---------------------- 
    157          zfact = 0.5e0 
     158         zfact = 0.5_wp 
    158159         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:) 
    159160      ENDIF 
    … …  
    226227      ENDIF 
    227228  
    228       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
     229      IF( l_trdtra )   THEN                      ! save the trends for further diagnostics 
    229230         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    230231         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    231          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_nsr, ztrdt ) 
    232          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_nsr, ztrds ) 
     232         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt ) 
     233         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds ) 
    233234         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )  
    234235      ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r3294 r3876  
    1919   USE sbc_oce         ! surface boundary condition: ocean 
    2020   USE dynspg_oce 
    21  
    2221   USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp     routine) 
    2322   USE trazdf_imp      ! vertical diffusion: implicit (tra_zdf_imp     routine) 
    24  
    2523   USE ldftra_oce      ! ocean active tracers: lateral physics 
    26    USE trdmod_oce      ! ocean active tracers: lateral physics 
    27    USE trdtra      ! ocean tracers trends  
     24   USE trd_oce         ! trends: ocean variables 
     25   USE trdtra          ! trends manager: tracers  
     26   ! 
    2827   USE in_out_manager  ! I/O manager 
    2928   USE prtctl          ! Print control 
    … …  
    9695            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk) 
    9796         END DO 
    98          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdf, ztrdt ) 
    99          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdf, ztrds ) 
     97         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
     98         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    10099         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    101100      ENDIF 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90

    r3632 r3876  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  trdtra  *** 
    4    !! Ocean diagnostics:  ocean tracers trends 
     4   !! Ocean diagnostics:  ocean tracers trends pre-processing 
    55   !!===================================================================== 
    6    !! History :  1.0  !  2004-08  (C. Talandier) Original code 
    7    !!            2.0  !  2005-04  (C. Deltel)    Add Asselin trend in the ML budget 
    8    !!            3.3  !  2010-06  (C. Ethe) merge TRA-TRC  
    9    !!---------------------------------------------------------------------- 
    10 #if  defined key_trdtra || defined key_trdtrc || defined key_trdmld || defined key_trdmld_trc  
    11    !!---------------------------------------------------------------------- 
    12    !!   trd_tra      : Call the trend to be computed 
    13    !!---------------------------------------------------------------------- 
    14    USE dom_oce          ! ocean domain  
    15    USE trdmod_oce       ! ocean active mixed layer tracers trends  
    16    USE trdmod           ! ocean active mixed layer tracers trends  
    17    USE trdmod_trc       ! ocean passive mixed layer tracers trends  
    18    USE in_out_manager   ! I/O manager 
    19    USE lib_mpp          ! MPP library 
    20    USE wrk_nemo        ! Memory allocation 
    21  
     6   !! History :  3.3  !  2010-06  (C. Ethe) creation for the TRA/TRC merge 
     7   !!            3.5  !  2012-02  (G. Madec) update the comments  
     8   !!---------------------------------------------------------------------- 
     9 
     10   !!---------------------------------------------------------------------- 
     11   !!   trd_tra       : pre-process the tracer trends 
     12   !!   trd_tra_adv   : transform a div(U.T) trend into a U.grad(T) trend 
     13   !!   trd_tra_mng   : tracer trend manager: dispatch to the diagnostic modules 
     14   !!   trd_tra_iom   : output 3D tracer trends using IOM 
     15   !!---------------------------------------------------------------------- 
     16   USE oce            ! ocean dynamics and tracers variables 
     17   USE dom_oce        ! ocean domain  
     18   USE sbc_oce        ! surface boundary condition: ocean 
     19   USE zdf_oce        ! ocean vertical physics 
     20   USE trd_oce        ! trends: ocean variables 
     21   USE trdmod_trc     ! ocean passive mixed layer tracers trends  
     22   USE trdglo         ! trends: global domain averaged 
     23   USE trdpen         ! trends: Potential ENergy 
     24   USE trdmxl         ! ocean active mixed layer tracers trends  
     25   USE ldftra_oce     ! ocean active tracers lateral physics 
     26   USE zdfddm         ! vertical physics: double diffusion 
     27   USE phycst         ! physical constants 
     28   USE in_out_manager ! I/O manager 
     29   USE iom            ! I/O manager library 
     30   USE lib_mpp        ! MPP library 
     31   USE wrk_nemo       ! Memory allocation 
    2232 
    2333   IMPLICIT NONE 
    2434   PRIVATE 
    2535 
    26    PUBLIC   trd_tra          ! called by all  traXX modules 
    27   
    28    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: trdtx, trdty, trdt  !: 
     36   PUBLIC   trd_tra   ! called by all tra_... modules 
     37 
     38   REAL(wp) ::   r2dt   ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0 
     39 
     40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   trdtx, trdty, trdt   ! use to store the temperature trends 
    2941 
    3042   !! * Substitutions 
    3143#  include "domzgr_substitute.h90" 
     44#  include "zdfddm_substitute.h90" 
    3245#  include "vectopt_loop_substitute.h90" 
    3346   !!---------------------------------------------------------------------- 
    34    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     47   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    3548   !! $Id$ 
    3649   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    3952 
    4053   INTEGER FUNCTION trd_tra_alloc() 
    41       !!---------------------------------------------------------------------------- 
     54      !!--------------------------------------------------------------------- 
    4255      !!                  ***  FUNCTION trd_tra_alloc  *** 
    43       !!---------------------------------------------------------------------------- 
     56      !!--------------------------------------------------------------------- 
    4457      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc ) 
    4558      ! 
    … …  
    5366      !!                  ***  ROUTINE trd_tra  *** 
    5467      !!  
    55       !! ** Purpose : Dispatch all trends computation, e.g. vorticity, mld or  
    56       !!              integral constraints 
     68      !! ** Purpose : pre-process tracer trends 
    5769      !! 
    58       !! ** Method/usage : For the mixed-layer trend, the control surface can be either 
    59       !!       a mixed layer depth (time varying) or a fixed surface (jk level or bowl).  
    60       !!      Choose control surface with nn_ctls in namelist NAMTRD : 
    61       !!        nn_ctls = 0  : use mixed layer with density criterion  
    62       !!        nn_ctls = 1  : read index from file 'ctlsurf_idx' 
    63       !!        nn_ctls > 1  : use fixed level surface jk = nn_ctls 
    64       !!---------------------------------------------------------------------- 
    65       ! 
    66       INTEGER                         , INTENT(in)           ::  kt      ! time step 
    67       CHARACTER(len=3)                , INTENT(in)           ::  ctype   ! tracers trends type 'TRA'/'TRC' 
    68       INTEGER                         , INTENT(in)           ::  ktra    ! tracer index 
    69       INTEGER                         , INTENT(in)           ::  ktrd    ! tracer trend index 
    70       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)           ::  ptrd    ! tracer trend  or flux 
    71       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::  pun     ! velocity  
    72       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::  ptra    ! Tracer variablea 
    73       ! 
    74       REAL(wp), POINTER, DIMENSION(:,:,:)  ::  ztrds 
    75       !!---------------------------------------------------------------------- 
    76  
     70      !! ** Method  : - mask the trend 
     71      !!              - advection (ptra present) converte the incoming flux (U.T)  
     72      !!              into trend (U.T => -U.grat(T)=div(U.T)-T.div(U)) through a  
     73      !!              call to trd_tra_adv 
     74      !!              - 'TRA' case : regroup T & S trends 
     75      !!              - send the trends to trd_tra_mng (trdmod_trc) for further processing 
     76      !!---------------------------------------------------------------------- 
     77      INTEGER                         , INTENT(in)           ::   kt      ! time step 
     78      CHARACTER(len=3)                , INTENT(in)           ::   ctype   ! tracers trends type 'TRA'/'TRC' 
     79      INTEGER                         , INTENT(in)           ::   ktra    ! tracer index 
     80      INTEGER                         , INTENT(in)           ::   ktrd    ! tracer trend index 
     81      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)           ::   ptrd    ! tracer trend  or flux 
     82      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pun     ! now velocity  
     83      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   ptra    ! now tracer variable 
     84      ! 
     85      INTEGER  ::   jk   ! loop indices 
     86      REAL(wp), POINTER, DIMENSION(:,:,:)  ::   zwt, zws, ztrdt, ztrds   ! 3D workspace 
     87      !!---------------------------------------------------------------------- 
     88      ! 
    7789      CALL wrk_alloc( jpi, jpj, jpk, ztrds ) 
    78  
    79       IF( .NOT. ALLOCATED( trdtx ) ) THEN       ! allocate trdtra arrays 
     90      !       
     91      IF( .NOT. ALLOCATED( trdtx ) ) THEN      ! allocate trdtra arrays 
    8092         IF( trd_tra_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' ) 
    8193      ENDIF 
    82        
    83       ! Control of optional arguments 
    84       IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN  
    85          IF( PRESENT( ptra ) ) THEN     
    86             SELECT CASE( ktrd )            ! shift depending on the direction 
    87             CASE( jptra_trd_xad )  ;  CALL trd_tra_adv( ptrd, pun, ptra, 'X', trdtx )  
    88             CASE( jptra_trd_yad )  ;  CALL trd_tra_adv( ptrd, pun, ptra, 'Y', trdty )  
    89             CASE( jptra_trd_zad )  ;  CALL trd_tra_adv( ptrd, pun, ptra, 'Z', trdt  )  
    90             END SELECT 
    91          ELSE 
    92             trdt(:,:,:) = ptrd(:,:,:) 
    93             IF( ktrd == jptra_trd_bbc .OR. ktrd == jptra_trd_qsr ) THEN 
    94                ztrds(:,:,:) = 0. 
    95                CALL trd_mod( trdt, ztrds, ktrd, ctype, kt ) 
    96             END IF 
    97          END IF 
    98       END IF 
    99  
    100       IF( ctype == 'TRA' .AND. ktra == jp_sal ) THEN  
    101          IF( PRESENT( ptra ) ) THEN     
    102             SELECT CASE( ktrd )            ! shift depending on the direction 
    103             CASE( jptra_trd_xad )   
    104                                 CALL trd_tra_adv( ptrd, pun, ptra, 'X', ztrds )  
    105                                 CALL trd_mod( trdtx, ztrds, ktrd, ctype, kt   ) 
    106             CASE( jptra_trd_yad )   
    107                                 CALL trd_tra_adv( ptrd, pun, ptra, 'Y', ztrds )  
    108                                 CALL trd_mod( trdty, ztrds, ktrd, ctype, kt   ) 
    109             CASE( jptra_trd_zad )     
    110                                 CALL trd_tra_adv( ptrd, pun, ptra, 'Z', ztrds )  
    111                                 CALL trd_mod( trdt , ztrds, ktrd, ctype, kt   ) 
    112             END SELECT 
    113          ELSE 
    114             ztrds(:,:,:) = ptrd(:,:,:) 
    115             CALL trd_mod( trdt, ztrds, ktrd, ctype, kt )   
    116          END IF 
    117       END IF 
    118  
    119       IF( ctype == 'TRC' ) THEN 
    120          ! 
    121          IF( PRESENT( ptra ) ) THEN   
    122             SELECT CASE( ktrd )            ! shift depending on the direction 
    123             CASE( jptra_trd_xad )   
    124                                 CALL trd_tra_adv( ptrd, pun, ptra, 'X', ztrds )  
    125                                 CALL trd_mod_trc( ztrds, ktra, ktrd, kt       ) 
    126             CASE( jptra_trd_yad )   
    127                                 CALL trd_tra_adv( ptrd, pun, ptra, 'Y', ztrds )  
    128                                 CALL trd_mod_trc( ztrds, ktra, ktrd, kt       ) 
    129             CASE( jptra_trd_zad )     
    130                                 CALL trd_tra_adv( ptrd, pun, ptra, 'Z', ztrds )  
    131                                 CALL trd_mod_trc( ztrds, ktra, ktrd, kt       ) 
    132             END SELECT 
    133          ELSE 
    134             ztrds(:,:,:) = ptrd(:,:,:) 
    135             CALL trd_mod_trc( ztrds, ktra, ktrd, kt       )   
    136          END IF 
     94 
     95      IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN   !==  Temperature trend  ==! 
     96         ! 
     97         SELECT CASE( ktrd ) 
     98         !                            ! advection: transform the advective flux into a trend 
     99         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'X', trdtx )  
     100         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'Y', trdty )  
     101         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'Z', trdt  )  
     102         CASE( jptra_bbc,    &        ! qsr, bbc: on temperature only, send to trd_tra_mng 
     103            &  jptra_qsr )   ;   trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) 
     104                                 ztrds(:,:,:) = 0._wp 
     105                                 CALL trd_tra_mng( trdt, ztrds, ktrd, kt ) 
     106         CASE DEFAULT                 ! other trends: masked trends 
     107            trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)              ! mask & store 
     108         END SELECT 
     109         ! 
     110      ENDIF 
     111 
     112      IF( ctype == 'TRA' .AND. ktra == jp_sal ) THEN      !==  Salinity trends  ==! 
     113         ! 
     114         SELECT CASE( ktrd ) 
     115         !                            ! advection: transform the advective flux into a trend 
     116         !                            !            and send T & S trends to trd_tra_mng 
     117         CASE( jptra_xad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'X'  , ztrds )  
     118                                  CALL trd_tra_mng( trdtx, ztrds, ktrd, kt   ) 
     119         CASE( jptra_yad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'Y'  , ztrds )  
     120                                  CALL trd_tra_mng( trdty, ztrds, ktrd, kt   ) 
     121         CASE( jptra_zad  )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'Z'  , ztrds )  
     122                                  CALL trd_tra_mng( trdt , ztrds, ktrd, kt   ) 
     123         CASE( jptra_zdfp )           ! diagnose the "PURE" Kz trend (here: just before the swap) 
     124            !                         ! iso-neutral diffusion case otherwise jptra_zdf is "PURE" 
     125            CALL wrk_alloc( jpi, jpj, jpk, zwt, zws, ztrdt ) 
     126            ! 
     127            zwt(:,:, 1 ) = 0._wp   ;   zws(:,:, 1 ) = 0._wp            ! vertical diffusive fluxes 
     128            zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp 
     129            DO jk = 2, jpk 
     130               zwt(:,:,jk) =   avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
     131               zws(:,:,jk) = fsavs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
     132            END DO 
     133            ! 
     134            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp 
     135            DO jk = 1, jpkm1 
     136               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk) 
     137               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / fse3t(:,:,jk)  
     138            END DO 
     139            CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )   
     140            ! 
     141            CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt ) 
     142            ! 
     143         CASE DEFAULT                 ! other trends: mask and send T & S trends to trd_tra_mng 
     144            ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) 
     145            CALL trd_tra_mng( trdt, ztrds, ktrd, kt )   
     146         END SELECT 
     147      ENDIF 
     148 
     149      IF( ctype == 'TRC' ) THEN                           !==  passive tracer trend  ==! 
     150         ! 
     151         SELECT CASE( ktrd ) 
     152         !                            ! advection: transform the advective flux into a masked trend 
     153         CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'X', ztrds )  
     154         CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'Y', ztrds )  
     155         CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'Z', ztrds )  
     156         CASE DEFAULT                 ! other trends: just masked  
     157                                 ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) 
     158         END SELECT 
     159         !                            ! send trend to trd_mod_trc 
     160         CALL trd_mod_trc( ztrds, ktra, ktrd, kt )  
    137161         ! 
    138162      ENDIF 
    … …  
    147171      !!                  ***  ROUTINE trd_tra_adv  *** 
    148172      !!  
    149       !! ** Purpose :   transformed the i-, j- or k-advective flux into thes 
    150       !!              i-, j- or k-advective trends, resp. 
    151       !! ** Method  :   i-advective trends = -un. di-1[T] = -( di-1[fi] - tn di-1[un] ) 
    152       !!                k-advective trends = -un. di-1[T] = -( dj-1[fi] - tn dj-1[un] ) 
    153       !!                k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] ) 
    154       !!---------------------------------------------------------------------- 
    155       REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk) ::   pf      ! advective flux in one direction 
    156       REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk) ::   pun     ! now velocity  in one direction 
    157       REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk) ::   ptn     ! now or before tracer  
    158       CHARACTER(len=1), INTENT(in )                         ::   cdir    ! X/Y/Z direction 
    159       REAL(wp)        , INTENT(out), DIMENSION(jpi,jpj,jpk) ::   ptrd    ! advective trend in one direction 
     173      !! ** Purpose :   transformed a advective flux into a masked advective trends 
     174      !! 
     175      !! ** Method  :   use the following transformation: -div(U.T) = - U grad(T) + T.div(U) 
     176      !!       i-advective trends = -un. di-1[T] = -( di-1[fi] - tn di-1[un] ) 
     177      !!       j-advective trends = -un. di-1[T] = -( dj-1[fi] - tn dj-1[un] ) 
     178      !!       k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] ) 
     179      !!                where fi is the incoming advective flux. 
     180      !!---------------------------------------------------------------------- 
     181      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pf      ! advective flux in one direction 
     182      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pun     ! now velocity   in one direction 
     183      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptn     ! now or before tracer  
     184      CHARACTER(len=1)                , INTENT(in   ) ::   cdir    ! X/Y/Z direction 
     185      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   ptrd    ! advective trend in one direction 
    160186      ! 
    161187      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    162       INTEGER  ::   ii, ij, ik   ! index shift function of the direction 
    163       REAL(wp) ::   zbtr         ! local scalar 
    164       !!---------------------------------------------------------------------- 
    165  
    166       SELECT CASE( cdir )            ! shift depending on the direction 
    167       CASE( 'X' )   ;   ii = 1   ; ij = 0   ;   ik = 0      ! i-advective trend 
    168       CASE( 'Y' )   ;   ii = 0   ; ij = 1   ;   ik = 0      ! j-advective trend 
    169       CASE( 'Z' )   ;   ii = 0   ; ij = 0   ;   ik =-1      ! k-advective trend 
     188      INTEGER  ::   ii, ij, ik   ! index shift as function of the direction 
     189      !!---------------------------------------------------------------------- 
     190      ! 
     191      SELECT CASE( cdir )      ! shift depending on the direction 
     192      CASE( 'X' )   ;   ii = 1   ;   ij = 0   ;   ik = 0      ! i-trend 
     193      CASE( 'Y' )   ;   ii = 0   ;   ij = 1   ;   ik = 0      ! j-trend 
     194      CASE( 'Z' )   ;   ii = 0   ;   ij = 0   ;   ik =-1      ! k-trend 
    170195      END SELECT 
    171  
    172       !                              ! set to zero uncomputed values 
    173       ptrd(jpi,:,:) = 0.e0   ;   ptrd(1,:,:) = 0.e0 
    174       ptrd(:,jpj,:) = 0.e0   ;   ptrd(:,1,:) = 0.e0 
    175       ptrd(:,:,jpk) = 0.e0 
    176       ! 
    177       ! 
    178       DO jk = 1, jpkm1 
     196      ! 
     197      !                        ! set to zero uncomputed values 
     198      ptrd(jpi,:,:) = 0._wp   ;   ptrd(1,:,:) = 0._wp 
     199      ptrd(:,jpj,:) = 0._wp   ;   ptrd(:,1,:) = 0._wp 
     200      ptrd(:,:,jpk) = 0._wp 
     201      ! 
     202      DO jk = 1, jpkm1         ! advective trend 
    179203         DO jj = 2, jpjm1 
    180204            DO ji = fs_2, fs_jpim1   ! vector opt. 
    181                zbtr    = 1.e0/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    182                ptrd(ji,jj,jk) = - zbtr * (      pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                    & 
    183                  &                          - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk)  ) 
     205               ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        & 
     206                 &                  - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk)  )   & 
     207                 &              / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )  * tmask(ji,jj,jk) 
    184208            END DO 
    185209         END DO 
    … …  
    188212   END SUBROUTINE trd_tra_adv 
    189213 
    190 #   else 
    191    !!---------------------------------------------------------------------- 
    192    !!   Default case :          Dummy module           No trend diagnostics 
    193    !!---------------------------------------------------------------------- 
    194    USE par_oce      ! ocean variables trends 
    195 CONTAINS 
    196    SUBROUTINE trd_tra( kt, ctype, ktra, ktrd, ptrd, pu, ptra ) 
    197       !!---------------------------------------------------------------------- 
    198       INTEGER                         , INTENT(in)           ::  kt      ! time step 
    199       CHARACTER(len=3)                , INTENT(in)           ::  ctype   ! tracers trends type 'TRA'/'TRC' 
    200       INTEGER                         , INTENT(in)           ::  ktra    ! tracer index 
    201       INTEGER                         , INTENT(in)           ::  ktrd    ! tracer trend index 
    202       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)           ::  ptrd    ! tracer trend  
    203       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::  pu      ! velocity  
    204       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::  ptra    ! Tracer variable  
    205       WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd(1,1,1), ptra(1,1,1), pu(1,1,1),   & 
    206          &                                                               ktrd, ktra, ctype, kt 
    207    END SUBROUTINE trd_tra 
    208 #   endif 
     214 
     215   SUBROUTINE trd_tra_mng( ptrdx, ptrdy, ktrd, kt ) 
     216      !!--------------------------------------------------------------------- 
     217      !!                  ***  ROUTINE trd_tra_mng  *** 
     218      !!  
     219      !! ** Purpose :   Dispatch all tracer trends computation, e.g. 3D output, 
     220      !!                integral constraints, potential energy, and/or  
     221      !!                mixed layer budget. 
     222      !!---------------------------------------------------------------------- 
     223      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend  
     224      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend 
     225      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index 
     226      INTEGER                   , INTENT(in   ) ::   kt      ! time step 
     227      !!---------------------------------------------------------------------- 
     228 
     229      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restart with Euler time stepping) 
     230      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog) 
     231      ENDIF 
     232 
     233      !                   ! 3D output of tracers trends using IOM interface 
     234      IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt ) 
     235 
     236      !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     237      IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt ) 
     238 
     239      !                   ! Potential ENergy trends 
     240      IF( ln_PE_trd  )   CALL trd_pen( ptrdx, ptrdy, ktrd, kt, r2dt ) 
     241 
     242      !                   ! Mixed layer trends for active tracers 
     243      IF( ln_tra_mxl )   THEN    
     244         !----------------------------------------------------------------------------------------------- 
     245         ! W.A.R.N.I.N.G : 
     246         ! jptra_ldf : called by traldf.F90 
     247         !                 at this stage we store: 
     248         !                  - the lateral geopotential diffusion (here, lateral = horizontal) 
     249         !                  - and the iso-neutral diffusion if activated  
     250         ! jptra_zdf : called by trazdf.F90 
     251         !                 * in case of iso-neutral diffusion we store the vertical diffusion component in the  
     252         !                   lateral trend including the K_z contrib, which will be removed later (see trd_mxl) 
     253         !----------------------------------------------------------------------------------------------- 
     254 
     255         SELECT CASE ( ktrd ) 
     256         CASE ( jptra_xad )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_xad, '3D' )   ! zonal    advection 
     257         CASE ( jptra_yad )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_yad, '3D' )   ! merid.   advection 
     258         CASE ( jptra_zad )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_zad, '3D' )   ! vertical advection 
     259         CASE ( jptra_ldf )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_ldf, '3D' )   ! lateral  diffusion 
     260         CASE ( jptra_bbl )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_bbl, '3D' )   ! bottom boundary layer 
     261         CASE ( jptra_zdf ) 
     262            IF( ln_traldf_iso ) THEN ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_ldf, '3D' )   ! lateral  diffusion (K_z) 
     263            ELSE                   ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_zdf, '3D' )   ! vertical diffusion (K_z) 
     264            ENDIF 
     265         CASE ( jptra_dmp )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_dmp, '3D' )   ! internal 3D restoring (tradmp) 
     266         CASE ( jptra_qsr )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_for, '3D' )   ! air-sea : penetrative sol radiat 
     267         CASE ( jptra_nsr )        ;   ptrdx(:,:,2:jpk) = 0._wp   ;   ptrdy(:,:,2:jpk) = 0._wp 
     268                                       CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_for, '2D' )   ! air-sea : non penetr sol radiation 
     269         CASE ( jptra_bbc )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_bbc, '3D' )   ! bottom bound cond (geoth flux) 
     270         CASE ( jptra_npc )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_npc, '3D' )   ! non penetr convect adjustment 
     271         CASE ( jptra_atf )        ;   CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_atf, '3D' )   ! asselin time filter (last trend) 
     272                                   ! 
     273                                       CALL trd_mxl( kt, r2dt )                             ! trends: Mixed-layer (output) 
     274         END SELECT 
     275         ! 
     276      ENDIF 
     277      ! 
     278   END SUBROUTINE trd_tra_mng 
     279 
     280 
     281   SUBROUTINE trd_tra_iom( ptrdx, ptrdy, ktrd, kt ) 
     282      !!--------------------------------------------------------------------- 
     283      !!                  ***  ROUTINE trd_tra_iom  *** 
     284      !!  
     285      !! ** Purpose :   output 3D tracer trends using IOM 
     286      !!---------------------------------------------------------------------- 
     287      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend  
     288      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend 
     289      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index 
     290      INTEGER                   , INTENT(in   ) ::   kt      ! time step 
     291      !! 
     292      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     293      INTEGER ::   ikbu, ikbv   ! local integers 
     294      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace  
     295      !!---------------------------------------------------------------------- 
     296      ! 
     297!!gm Rq: mask the trends already masked in trd_tra, but lbc_lnk should probably be added 
     298      ! 
     299      SELECT CASE( ktrd ) 
     300      CASE( jptra_xad  )   ;   CALL iom_put( "ttrd_xad" , ptrdx )        ! x- horizontal advection 
     301                               CALL iom_put( "strd_xad" , ptrdy ) 
     302      CASE( jptra_yad  )   ;   CALL iom_put( "ttrd_yad" , ptrdx )        ! y- horizontal advection 
     303                               CALL iom_put( "strd_yad" , ptrdy ) 
     304      CASE( jptra_zad  )   ;   CALL iom_put( "ttrd_zad" , ptrdx )        ! z- vertical   advection 
     305                               CALL iom_put( "strd_zad" , ptrdy ) 
     306                               IF( .NOT. lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface 
     307                                  CALL wrk_alloc( jpi, jpj, z2dx, z2dy ) 
     308                                  z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1) 
     309                                  z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1) 
     310                                  CALL iom_put( "ttrd_sad", z2dx ) 
     311                                  CALL iom_put( "strd_sad", z2dy ) 
     312                                  CALL wrk_dealloc( jpi, jpj, z2dx, z2dy ) 
     313                               ENDIF 
     314      CASE( jptra_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion 
     315                               CALL iom_put( "strd_ldf" , ptrdy ) 
     316      CASE( jptra_zdf  )   ;   CALL iom_put( "ttrd_zdf" , ptrdx )        ! vertical diffusion (including Kz contribution) 
     317                               CALL iom_put( "strd_zdf" , ptrdy ) 
     318      CASE( jptra_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution) 
     319                               CALL iom_put( "strd_zdfp", ptrdy ) 
     320      CASE( jptra_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping) 
     321                               CALL iom_put( "strd_dmp" , ptrdy ) 
     322      CASE( jptra_bbl  )   ;   CALL iom_put( "ttrd_bbl" , ptrdx )        ! bottom boundary layer 
     323                               CALL iom_put( "strd_bbl" , ptrdy ) 
     324      CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing 
     325                               CALL iom_put( "strd_npc" , ptrdy ) 
     326      CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx )        ! surface forcing + runoff (ln_rnf=T) 
     327                               CALL iom_put( "strd_cdt" , ptrdy ) 
     328      CASE( jptra_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature) 
     329      CASE( jptra_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature) 
     330      CASE( jptra_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter 
     331                               CALL iom_put( "strd_atf" , ptrdy ) 
     332      END SELECT 
     333      ! 
     334   END SUBROUTINE trd_tra_iom 
     335 
    209336   !!====================================================================== 
    210337END MODULE trdtra 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90

    r3294 r3876  
    44   !! Ocean diagnostics:  momentum trends 
    55   !!===================================================================== 
    6    !! History :  1.0  !  04-2006  (L. Brunier, A-M. Treguier) Original code  
    7    !!            2.0  !  04-2008  (C. Talandier) New trends organization 
     6   !! History :  1.0  !  2006-01  (L. Brunier, A-M. Treguier) Original code  
     7   !!            2.0  !  2008-04  (C. Talandier) New trends organization 
     8   !!            3.5  !  2012-02  (G. Madec) regroup beta.V computation with pvo trend 
    89   !!---------------------------------------------------------------------- 
    9 #if defined key_trdvor   ||   defined key_esopa 
    10    !!---------------------------------------------------------------------- 
    11    !!   'key_trdvor'   : momentum trend diagnostics 
     10 
    1211   !!---------------------------------------------------------------------- 
    1312   !!   trd_vor      : momentum trends averaged over the depth 
    … …  
    1716   USE oce             ! ocean dynamics and tracers variables 
    1817   USE dom_oce         ! ocean space and time domain variables 
    19    USE trdmod_oce      ! ocean variables trends 
     18   USE trd_oce         ! trends: ocean variables 
    2019   USE zdf_oce         ! ocean vertical physics 
    21    USE in_out_manager  ! I/O manager 
     20   USE sbc_oce         ! surface boundary condition: ocean 
    2221   USE phycst          ! Define parameters for the routines 
    2322   USE ldfdyn_oce      ! ocean active tracers: lateral physics 
    2423   USE dianam          ! build the name of file (routine) 
    2524   USE zdfmxl          ! mixed layer depth 
     25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     26   USE in_out_manager  ! I/O manager 
    2627   USE ioipsl          ! NetCDF library 
    27    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2828   USE lib_mpp         ! MPP library 
    2929   USE wrk_nemo        ! Memory allocation 
    30  
    3130 
    3231   IMPLICIT NONE 
    … …  
    3736   END INTERFACE 
    3837 
    39    PUBLIC   trd_vor        ! routine called by step.F90 
    40    PUBLIC   trd_vor_zint   ! routine called by dynamics routines 
     38   PUBLIC   trd_vor        ! routine called by trddyn.F90 
    4139   PUBLIC   trd_vor_init   ! routine called by opa.F90 
    4240   PUBLIC   trd_vor_alloc  ! routine called by nemogcm.F90 
    … …  
    8078      IF( trd_vor_alloc /= 0 )   CALL ctl_warn('trd_vor_alloc: failed to allocate arrays') 
    8179   END FUNCTION trd_vor_alloc 
     80 
     81 
     82   SUBROUTINE trd_vor( putrd, pvtrd, ktrd, kt ) 
     83      !!---------------------------------------------------------------------- 
     84      !!                  ***  ROUTINE trd_vor  *** 
     85      !!  
     86      !! ** Purpose :  computation of cumulated trends over analysis period 
     87      !!               and make outputs (NetCDF or DIMG format) 
     88      !!---------------------------------------------------------------------- 
     89      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends  
     90      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index 
     91      INTEGER                   , INTENT(in   ) ::   kt             ! time step 
     92      ! 
     93      INTEGER ::   ji, jj   ! dummy loop indices 
     94      REAL(wp), POINTER, DIMENSION(:,:) ::   ztswu, ztswv    ! 2D workspace  
     95      !!---------------------------------------------------------------------- 
     96 
     97      CALL wrk_alloc( jpi, jpj, ztswu, ztswv ) 
     98 
     99      SELECT CASE( ktrd )  
     100      CASE( jpdyn_hpg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_prg )   ! Hydrostatique Pressure Gradient  
     101      CASE( jpdyn_keg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_keg )   ! KE Gradient  
     102      CASE( jpdyn_rvo )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_rvo )   ! Relative Vorticity  
     103      CASE( jpdyn_pvo )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_pvo )   ! Planetary Vorticity Term  
     104      CASE( jpdyn_ldf )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_ldf )   ! Horizontal Diffusion  
     105      CASE( jpdyn_zad )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_zad )   ! Vertical Advection  
     106      CASE( jpdyn_spg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_spg )   ! Surface Pressure Grad.  
     107      CASE( jpdyn_zdf )                                                      ! Vertical Diffusion  
     108         ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0 
     109         DO jj = 2, jpjm1                                                             ! wind stress trends 
     110            DO ji = fs_2, fs_jpim1   ! vector opt. 
     111               ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 ) 
     112               ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 ) 
     113            END DO 
     114         END DO 
     115         ! 
     116         CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf )                             ! zdf trend including surf./bot. stresses  
     117         CALL trd_vor_zint( ztswu, ztswv, jpvor_swf )                             ! surface wind stress  
     118      CASE( jpdyn_bfr ) 
     119         CALL trd_vor_zint( putrd, pvtrd, jpvor_bfr )                             ! Bottom stress 
     120         ! 
     121      CASE( jpdyn_atf )       ! last trends: perform the output of 2D vorticity trends 
     122         CALL trd_vor_iom( kt ) 
     123      END SELECT 
     124      ! 
     125      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv ) 
     126      ! 
     127   END SUBROUTINE trd_vor 
    82128 
    83129 
    … …  
    109155      !!      trends output in netCDF format using ioipsl 
    110156      !!---------------------------------------------------------------------- 
    111       ! 
    112157      INTEGER                     , INTENT(in   ) ::   ktrd       ! ocean trend index 
    113158      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   putrdvor   ! u vorticity trend  
    … …  
    131176      !  ===================================== 
    132177 
    133       SELECT CASE (ktrd)  
    134       ! 
    135       CASE (jpvor_bfr)        ! bottom friction 
     178      SELECT CASE( ktrd )  
     179      ! 
     180      CASE( jpvor_bfr )        ! bottom friction 
    136181         DO jj = 2, jpjm1 
    137182            DO ji = fs_2, fs_jpim1  
    … …  
    143188         END DO 
    144189         ! 
    145       CASE (jpvor_swf)        ! wind stress 
     190      CASE( jpvor_swf )        ! wind stress 
    146191         zudpvor(:,:) = putrdvor(:,:) * fse3u(:,:,1) * e1u(:,:) * umask(:,:,1) 
    147192         zvdpvor(:,:) = pvtrdvor(:,:) * fse3v(:,:,1) * e2v(:,:) * vmask(:,:,1) 
    … …  
    154199    
    155200      ! Curl 
    156       DO ji=1,jpim1 
    157          DO jj=1,jpjm1 
     201      DO ji = 1, jpim1 
     202         DO jj = 1, jpjm1 
    158203            vortrd(ji,jj,ktrd) = (    zvdpvor(ji+1,jj) - zvdpvor(ji,jj)       & 
    159204                 &                - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) )   ) / ( e1f(ji,jj) * e2f(ji,jj) ) 
    … …  
    229274      END DO 
    230275 
    231       ! Save Beta.V term to avoid average before Curl 
    232       ! Beta.V : intergration, no average 
    233       IF( ktrd == jpvor_bev ) THEN  
     276      ! Planetary vorticity: 2nd computation (Beta.V term) store the vertical sum 
     277      ! as Beta.V term need intergration, not average 
     278      IF( ktrd == jpvor_pvo ) THEN  
    234279         zubet(:,:) = zudpvor(:,:) 
    235280         zvbet(:,:) = zvdpvor(:,:) 
    236       ENDIF 
    237  
    238       ! Average except for Beta.V 
     281         DO ji = 1, jpim1 
     282            DO jj = 1, jpjm1 
     283               vortrd(ji,jj,jpvor_bev) = (    zvbet(ji+1,jj) - zvbet(ji,jj)     & 
     284                  &                       - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) 
     285            END DO 
     286         END DO 
     287         ! Average of the Curl and Surface mask 
     288         vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * hur(:,:) * fmask(:,:,1) 
     289      ENDIF 
     290      ! 
     291      ! Average  
    239292      zudpvor(:,:) = zudpvor(:,:) * hur(:,:) 
    240293      zvdpvor(:,:) = zvdpvor(:,:) * hvr(:,:) 
    241     
     294      ! 
    242295      ! Curl 
    243296      DO ji=1,jpim1 
    … …  
    247300         END DO 
    248301      END DO 
    249  
    250302      ! Surface mask 
    251303      vortrd(:,:,ktrd) = vortrd(:,:,ktrd) * fmask(:,:,1) 
    252  
    253       ! Special treatement for the Beta.V term 
    254       ! Compute the Curl of the Beta.V term which is not averaged 
    255       IF( ktrd == jpvor_bev ) THEN 
    256          DO ji=1,jpim1 
    257             DO jj=1,jpjm1 
    258                vortrd(ji,jj,jpvor_bev) = (    zvbet(ji+1,jj) - zvbet(ji,jj)     & 
    259                   &                       - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) 
    260             END DO 
    261          END DO 
    262  
    263          ! Average on the Curl 
    264          vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * hur(:,:) 
    265  
    266          ! Surface mask 
    267          vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * fmask(:,:,1) 
    268       ENDIF 
    269304    
    270305      IF( ndebug /= 0 ) THEN 
    … …  
    278313 
    279314 
    280    SUBROUTINE trd_vor( kt ) 
     315   SUBROUTINE trd_vor_iom( kt ) 
    281316      !!---------------------------------------------------------------------- 
    282317      !!                  ***  ROUTINE trd_vor  *** 
    … …  
    285320      !!               and make outputs (NetCDF or DIMG format) 
    286321      !!---------------------------------------------------------------------- 
    287       ! 
    288       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     322      INTEGER                   , INTENT(in   ) ::   kt             ! time step 
    289323      ! 
    290324      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    … …  
    305339 
    306340      IF( kt > nit000 )   vor_avrb(:,:) = vor_avr(:,:) 
    307  
    308       IF( ndebug /= 0 ) THEN 
    309           WRITE(numout,*) ' debuging trd_vor: I.1 done ' 
    310           CALL FLUSH(numout) 
    311       ENDIF 
    312341 
    313342      ! I.2 vertically integrated vorticity 
    … …  
    330359 
    331360      ! Curl 
    332       DO ji=1,jpim1 
    333          DO jj=1,jpjm1 
     361      DO ji = 1, jpim1 
     362         DO jj = 1, jpjm1 
    334363            vor_avr(ji,jj) = (  ( zvn(ji+1,jj) - zvn(ji,jj) )    & 
    335364               &              - ( zun(ji,jj+1) - zun(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1) 
    … …  
    337366      END DO 
    338367       
    339       IF( ndebug /= 0 ) THEN 
    340          WRITE(numout,*) ' debuging trd_vor: I.2 done' 
    341          CALL FLUSH(numout) 
    342       ENDIF 
    343  
    344368      !  ================================= 
    345369      !   II. Cumulated trends 
    … …  
    351375         vor_avrbb(:,:) = vor_avrb(:,:) 
    352376         vor_avrbn(:,:) = vor_avr (:,:) 
    353       ENDIF 
    354  
    355       IF( ndebug /= 0 ) THEN 
    356          WRITE(numout,*) ' debuging trd_vor: I1.1 done' 
    357          CALL FLUSH(numout) 
    358377      ENDIF 
    359378 
    … …  
    371390      ENDIF 
    372391 
    373       IF( ndebug /= 0 ) THEN 
    374          WRITE(numout,*) ' debuging trd_vor: II.2 done' 
    375          CALL FLUSH(numout) 
    376       ENDIF 
    377  
    378392      !  ============================================= 
    379393      !   III. Output in netCDF + residual computation 
    … …  
    391405         vor_avrtot(:,:) = (  vor_avr(:,:) - vor_avrbn(:,:) + vor_avrb(:,:) - vor_avrbb(:,:) ) * zmean 
    392406 
    393          IF( ndebug /= 0 ) THEN 
    394              WRITE(numout,*) ' zmean = ',zmean 
    395              WRITE(numout,*) ' debuging trd_vor: III.1 done' 
    396              CALL FLUSH(numout) 
    397          ENDIF 
    398407 
    399408         ! III.2 compute residual 
    … …  
    406415         CALL lbc_lnk( vor_avrres, 'F', 1. ) 
    407416 
    408          IF( ndebug /= 0 ) THEN 
    409             WRITE(numout,*) ' debuging trd_vor: III.2 done' 
    410             CALL FLUSH(numout) 
    411          ENDIF 
    412417 
    413418         ! III.3 time evolution array swap 
    … …  
    415420         vor_avrbb(:,:) = vor_avrb(:,:) 
    416421         vor_avrbn(:,:) = vor_avr (:,:) 
    417  
    418          IF( ndebug /= 0 ) THEN 
    419             WRITE(numout,*) ' debuging trd_vor: III.3 done' 
    420             CALL FLUSH(numout) 
    421          ENDIF 
    422422         ! 
    423423         nmoydpvor = 0 
    … …  
    463463      CALL wrk_dealloc( jpi, jpj, zun, zvn )                                    
    464464      ! 
    465    END SUBROUTINE trd_vor 
     465   END SUBROUTINE trd_vor_iom 
    466466 
    467467 
    … …  
    587587   END SUBROUTINE trd_vor_init 
    588588 
    589 #else 
    590    !!---------------------------------------------------------------------- 
    591    !!   Default option :                                       Empty module 
    592    !!---------------------------------------------------------------------- 
    593    INTERFACE trd_vor_zint 
    594       MODULE PROCEDURE trd_vor_zint_2d, trd_vor_zint_3d 
    595    END INTERFACE 
    596 CONTAINS 
    597    SUBROUTINE trd_vor( kt )        ! Empty routine 
    598       WRITE(*,*) 'trd_vor: You should not have seen this print! error?', kt 
    599    END SUBROUTINE trd_vor 
    600    SUBROUTINE trd_vor_zint_2d( putrdvor, pvtrdvor, ktrd ) 
    601       REAL, DIMENSION(:,:), INTENT( inout ) ::   putrdvor, pvtrdvor 
    602       INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index 
    603       WRITE(*,*) 'trd_vor_zint_2d: You should not have seen this print! error?', putrdvor(1,1), pvtrdvor(1,1), ktrd 
    604    END SUBROUTINE trd_vor_zint_2d 
    605    SUBROUTINE trd_vor_zint_3d( putrdvor, pvtrdvor, ktrd ) 
    606       REAL, DIMENSION(:,:,:), INTENT( inout ) ::   putrdvor, pvtrdvor 
    607       INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index 
    608       WRITE(*,*) 'trd_vor_zint_3d: You should not have seen this print! error?', putrdvor(1,1,1), pvtrdvor(1,1,1), ktrd 
    609    END SUBROUTINE trd_vor_zint_3d 
    610    SUBROUTINE trd_vor_init              ! Empty routine 
    611       WRITE(*,*) 'trd_vor_init: You should not have seen this print! error?' 
    612    END SUBROUTINE trd_vor_init 
    613 #endif 
    614589   !!====================================================================== 
    615590END MODULE trdvor 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor_oce.F90

    r2715 r3876  
    44   !! Ocean trends :   set vorticity trend variables 
    55   !!====================================================================== 
    6    !! History :  9.0  !  04-2006  (L. Brunier, A-M. Treguier) Original code  
     6   !! History :  1.0  !  04-2006  (L. Brunier, A-M. Treguier) Original code  
    77   !!---------------------------------------------------------------------- 
    8  
    9    !!---------------------------------------------------------------------- 
     8    
    109   USE par_oce      ! ocean parameters 
    1110 
    … …  
    1312   PRIVATE 
    1413 
    15 #if defined key_trdvor 
    16    LOGICAL, PUBLIC, PARAMETER ::   lk_trdvor = .TRUE.    !: momentum trend flag 
    17 #else 
    18    LOGICAL, PUBLIC, PARAMETER ::   lk_trdvor = .FALSE.   !: momentum trend flag 
    19 #endif 
    2014   !                                               !!* vorticity trends index 
    2115   INTEGER, PUBLIC, PARAMETER ::   jpltot_vor = 11  !: Number of vorticity trend terms 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r3294 r3876  
    66   !! History :  1.0  ! 2003-08  (G. Madec)  original code 
    77   !!            3.2  ! 2009-07  (S. Masson, G. Madec)  IOM + merge of DO-loop 
     8   !!            3.5  ! 2012-03  (G. Madec)  make public the density criteria for trdmxl  
    89   !!---------------------------------------------------------------------- 
    910   !!   zdf_mxl      : Compute the turbocline and mixed layer depths. 
    … …  
    2930   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3031   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     32 
     33   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
     34   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
    3135 
    3236   !! * Substitutions 
    … …  
    6367      !! ** Method  :   The mixed layer depth is the shallowest W depth with  
    6468      !!      the density of the corresponding T point (just bellow) bellow a 
    65       !!      given value defined locally as rho(10m) + zrho_c 
     69      !!      given value defined locally as rho(10m) + rho_c 
    6670      !!               The turbocline depth is the depth at which the vertical 
    6771      !!      eddy diffusivity coefficient (resulting from the vertical physics 
    6872      !!      alone, not the isopycnal part, see trazdf.F) fall below a given 
    69       !!      value defined locally (avt_c here taken equal to 5 cm/s2) 
     73      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default) 
    7074      !! 
    7175      !! ** Action  :   nmln, hmld, hmlp, hmlpt 
    … …  
    7680      INTEGER  ::   iikn, iiki          ! temporary integer within a do loop 
    7781      INTEGER, POINTER, DIMENSION(:,:) ::   imld                ! temporary workspace 
    78       REAL(wp) ::   zrho_c = 0.01_wp    ! density criterion for mixed layer depth 
    79       REAL(wp) ::   zavt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
    8082      !!---------------------------------------------------------------------- 
    8183      ! 
    … …  
    98100         DO jj = 1, jpj 
    99101            DO ji = 1, jpi 
    100                IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + zrho_c )   nmln(ji,jj) = jk      ! Mixed layer 
    101                IF( avt (ji,jj,jk) < zavt_c                     )   imld(ji,jj) = jk      ! Turbocline  
     102               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rho_c )   nmln(ji,jj) = jk      ! Mixed layer 
     103               IF( avt (ji,jj,jk) < avt_c                     )   imld(ji,jj) = jk      ! Turbocline  
    102104            END DO 
    103105         END DO 

  • branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r3769 r3876  
    5959   USE zdfini          ! vertical physics setting          (zdf_init routine) 
    6060   USE phycst          ! physical constant                  (par_cst routine) 
    61    USE trdmod          ! momentum/tracers trends       (trd_mod_init routine) 
     61   USE trdini          ! dyn/tra trends initialization     (trd_init routine) 
    6262   USE asminc          ! assimilation increments      
    6363   USE asmbkg          ! writing out state trajectory 
    … …  
    395395      IF( lk_diadct     )   CALL dia_dct_init   ! Sections tranports 
    396396                            CALL dia_hsb_init   ! heat content, salt content and volume budgets 
    397                             CALL trd_mod_init   ! Mixed-layer/Vorticity/Integral constraints trends 
     397                            CALL     trd_init   ! momentum/tracer trends 
    398398      IF( lk_diaobs     ) THEN                  ! Observation & model comparison 
    399399                            CALL dia_obs_init            ! Initialize observational data 
    … …  
    574574      !! ** Method  : 
    575575      !!---------------------------------------------------------------------- 
    576       INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have 
     576      INTEGER, INTENT(in) ::   num_pes  ! The number of MPI processes we have 
    577577      ! 
    578578      INTEGER, PARAMETER :: nfactmax = 20 
    … …  
    583583      INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors 
    584584      !!---------------------------------------------------------------------- 
    585  
     585      ! 
    586586      ierr = 0 
    587  
     587      ! 
    588588      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr ) 
    589  
     589      ! 
    590590      IF( nfact <= 1 ) THEN 
    591591         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed' 
    … …  
    629629      INTEGER, PARAMETER :: ntest = 14 
    630630      INTEGER :: ilfax(ntest) 
    631  
     631      ! 
    632632      ! lfax contains the set of allowed factors. 
    633633      data (ilfax(jl),jl=1,ntest) / 16384, 8192, 4096, 2048, 1024, 512, 256,  & 
    … …  
    680680 
    681681#if defined key_mpp_mpi 
     682 
    682683   SUBROUTINE nemo_northcomms 
    683684      !!====================================================================== 
    684685      !!                     ***  ROUTINE  nemo_northcomms  *** 
    685       !! nemo_northcomms    :  Setup for north fold exchanges with explicit peer to peer messaging 
     686      !! nemo_northcomms :  Setup for north fold exchanges with explicit peer to peer messaging 
    686687      !!===================================================================== 
    687688      !!---------------------------------------------------------------------- 
    … …  
    691692      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE) 
    692693      !!---------------------------------------------------------------------- 
    693  
    694694      INTEGER ::   ji, jj, jk, ij, jtyp    ! dummy loop indices 
    695695      INTEGER ::   ijpj                    ! number of rows involved in north-fold exchange 
    696696      INTEGER ::   northcomms_alloc        ! allocate return status 
    697       REAL(wp), ALLOCATABLE, DIMENSION ( :,: ) ::   znnbrs     ! workspace 
    698       LOGICAL,  ALLOCATABLE, DIMENSION ( : )   ::   lrankset   ! workspace 
     697      REAL(wp), ALLOCATABLE, DIMENSION (:,:) ::   znnbrs     ! workspace 
     698      LOGICAL,  ALLOCATABLE, DIMENSION (:)   ::   lrankset   ! workspace 
     699      !!---------------------------------------------------------------------- 
    699700 
    700701      IF(lwp) WRITE(numout,*) 
    … …  
    702703      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    703704 
    704       !!---------------------------------------------------------------------- 
    705705      ALLOCATE( znnbrs(jpi,jpj), stat = northcomms_alloc ) 
    706706      ALLOCATE( lrankset(jpnij), stat = northcomms_alloc ) 
    … …  
    724724      ! Exchange and store ranks on northern rows 
    725725 
    726       DO jtyp = 1,4 
     726      DO jtyp = 1, 4 
    727727 
    728728         lrankset = .FALSE. 
    729          znnbrs = narea 
     729         znnbrs   = narea 
    730730         SELECT CASE (jtyp) 
    731731            CASE(1) 
    … …  
    739739         END SELECT 
    740740 
    741          IF ( njmppt(narea) .EQ. MAXVAL( njmppt ) ) THEN 
     741         IF ( njmppt(narea) == MAXVAL( njmppt ) ) THEN 
    742742            DO jj = nlcj-ijpj+1, nlcj 
    743743               ij = jj - nlcj + ijpj 
    … …  
    817817 
    818818   END SUBROUTINE nemo_northcomms 
     819    
    819820#else 
    820821   SUBROUTINE nemo_northcomms      ! Dummy routine 
    … …  
    822823   END SUBROUTINE nemo_northcomms 
    823824#endif 
     825 
    824826   !!====================================================================== 
    825827END MODULE nemogcm 
Note: See TracChangeset for help on using the changeset viewer.