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.
zdfddm.F90 in NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfddm.F90 @ 13228

Last change on this file since 13228 was 13228, checked in by smasson, 4 years ago

better e3: update with trunk@13227 see #2385
  • Property svn:keywords set to Id
File size: 8.7 KB
RevLine 
[3]1MODULE zdfddm
2   !!======================================================================
3   !!                       ***  MODULE  zdfddm  ***
4   !! Ocean physics : double diffusion mixing parameterization
5   !!======================================================================
[1601]6   !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
[4990]8   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
9   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta
[9019]10   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
[1601]11   !!----------------------------------------------------------------------
[9019]12
[3]13   !!----------------------------------------------------------------------
[9019]14   !!   zdf_ddm       : compute the Kz for salinity
[3]15   !!----------------------------------------------------------------------
[9019]16   USE oce            ! ocean dynamics and tracers variables
17   USE dom_oce        ! ocean space and time domain variables
18   USE zdf_oce        ! ocean vertical physics variables
[4990]19   USE eosbn2         ! equation of state
20   !
[9019]21   USE in_out_manager ! I/O manager
22   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
23   USE prtctl         ! Print control
24   USE lib_mpp        ! MPP library
[3]25
26   IMPLICIT NONE
27   PRIVATE
28
[2528]29   PUBLIC   zdf_ddm       ! called by step.F90
[3]30
31   !! * Substitutions
[12377]32#  include "do_loop_substitute.h90"
[3]33   !!----------------------------------------------------------------------
[9598]34   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]35   !! $Id$
[10068]36   !! Software governed by the CeCILL license (see ./LICENSE)
[3]37   !!----------------------------------------------------------------------
38CONTAINS
39
[12377]40   SUBROUTINE zdf_ddm( kt, Kmm, p_avm, p_avt, p_avs )
[2715]41      !!----------------------------------------------------------------------
[3]42      !!                  ***  ROUTINE zdf_ddm  ***
43      !!                   
44      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
[1601]45      !!              effect of salt fingering and diffusive convection.
[3]46      !!
47      !! ** Method  :   Diapycnal mixing is increased in case of double
48      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
49      !!      following Merryfield et al. (1999). The rate of double diffusive
[4990]50      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
[3]51      !!         * salt fingering (Schmitt 1981):
[4990]52      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
53      !!      for R > 1 and rn2 > 0 : zavfs = O
54      !!      otherwise                : zavft = 0.7 zavs / R
[3]55      !!         * diffusive layering (Federov 1988):
[4990]56      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
[3]57      !!      otherwise                   : zavdt = 0
[4990]58      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
59      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
[3]60      !!      otherwise                     : zavds = 0
61      !!         * update the eddy diffusivity:
62      !!      avt = avt + zavft + zavdt
63      !!      avs = avs + zavfs + zavds
[9019]64      !!      avm is required to remain at least above avt and avs.
[3]65      !!     
[1601]66      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
[3]67      !!
[1601]68      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
[3]69      !!----------------------------------------------------------------------
[12377]70      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
71      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
[9019]72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
74      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
[2715]75      !
[1601]76      INTEGER  ::   ji, jj , jk     ! dummy loop indices
[4990]77      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
78      REAL(wp) ::   zdt, zds
[13228]79      REAL(wp) ::   zinr            !   -      -
80      REAL(dp) ::         zrr       !   -      -
81      REAL(wp) ::   zavft           !   -      -
82      REAL(dp) ::          zavfs    !   -      -
[4990]83      REAL(wp) ::   zavdt, zavds    !   -      -
[9019]84      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
[3]85      !!----------------------------------------------------------------------
[3294]86      !
[3]87      !                                                ! ===============
88      DO jk = 2, jpkm1                                 ! Horizontal slab
89         !                                             ! ===============
90         ! Define the mask
91         ! ---------------
[9019]92!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
93!!gm                     ==>>>  test in the loop instead of use of mask arrays
94!!gm                            and many acces in memory
95         
[12377]96         DO_2D_11_11
97            zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
[12622]98!!gm please, use e3w at Kmm below
[12377]99               &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
100            !
101            zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
102                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
103            zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
104                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
105            !
106            zdt = zaw * ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )
107            zds = zbw * ( ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) 
108            IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
109            zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
110         END_2D
[3]111
[12377]112         DO_2D_11_11
113            ! stability indicator: msks=1 if rn2>0; 0 elsewhere
114            IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
115            ELSE                                       ;   zmsks(ji,jj) = 1._wp
116            ENDIF
117            ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
118            IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
119            ELSE                                       ;   zmskf(ji,jj) = 1._wp
120            ENDIF
121            ! diffusive layering indicators:
122            !     ! mskdl1=1 if 0< R <1; 0 elsewhere
123            IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
124            ELSE                                       ;   zmskd1(ji,jj) = 1._wp
125            ENDIF
126            !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
127            IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
128            ELSE                                       ;   zmskd2(ji,jj) = 1._wp
129            ENDIF
130            !   mskdl3=1 if 0.5< R <1; 0 elsewhere
131            IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
132            ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
133            ENDIF
134         END_2D
[3]135         ! mask zmsk in order to have avt and avs masked
[7753]136         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
[3]137
138
139         ! Update avt and avs
140         ! ------------------
141         ! Constant eddy coefficient: reset to the background value
[12377]142         DO_2D_11_11
143            zinr = 1._wp / zrau(ji,jj)
144            ! salt fingering
145            zrr = zrau(ji,jj) / rn_hsbfr
146            zrr = zrr * zrr
147            zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
148            zavft = 0.7 * zavfs * zinr
149            ! diffusive layering
150            zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
151            zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
152               &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
153            ! add to the eddy viscosity coef. previously computed
154            p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
155            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
156            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
157         END_2D
[3]158         !                                                ! ===============
159      END DO                                              !   End of slab
160      !                                                   ! ===============
[1601]161      !
[12377]162      IF(sn_cfctl%l_prtctl) THEN
[9440]163         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
[49]164      ENDIF
[1601]165      !
[3]166   END SUBROUTINE zdf_ddm
167   
168   !!======================================================================
169END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.