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_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfddm.F90 @ 13998

Last change on this file since 13998 was 13998, checked in by techene, 3 years ago

branch updated with trunk 13787
  • Property svn:keywords set to Id
File size: 8.8 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"
[13427]33#  include "domzgr_substitute.h90"
[3]34   !!----------------------------------------------------------------------
[9598]35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]36   !! $Id$
[10068]37   !! Software governed by the CeCILL license (see ./LICENSE)
[3]38   !!----------------------------------------------------------------------
39CONTAINS
40
[12377]41   SUBROUTINE zdf_ddm( kt, Kmm, p_avm, p_avt, p_avs )
[2715]42      !!----------------------------------------------------------------------
[3]43      !!                  ***  ROUTINE zdf_ddm  ***
44      !!                   
45      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
[1601]46      !!              effect of salt fingering and diffusive convection.
[3]47      !!
48      !! ** Method  :   Diapycnal mixing is increased in case of double
49      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
50      !!      following Merryfield et al. (1999). The rate of double diffusive
[4990]51      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
[3]52      !!         * salt fingering (Schmitt 1981):
[4990]53      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
54      !!      for R > 1 and rn2 > 0 : zavfs = O
55      !!      otherwise                : zavft = 0.7 zavs / R
[3]56      !!         * diffusive layering (Federov 1988):
[4990]57      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
[3]58      !!      otherwise                   : zavdt = 0
[4990]59      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
60      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
[3]61      !!      otherwise                     : zavds = 0
62      !!         * update the eddy diffusivity:
63      !!      avt = avt + zavft + zavdt
64      !!      avs = avs + zavfs + zavds
[9019]65      !!      avm is required to remain at least above avt and avs.
[3]66      !!     
[1601]67      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
[3]68      !!
[1601]69      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
[3]70      !!----------------------------------------------------------------------
[12377]71      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
72      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
[9019]73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
74      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
75      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
[2715]76      !
[1601]77      INTEGER  ::   ji, jj , jk     ! dummy loop indices
[4990]78      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
79      REAL(wp) ::   zdt, zds
[13226]80      REAL(wp) ::   zinr            !   -      -
81      REAL(dp) ::         zrr       !   -      -
82      REAL(wp) ::   zavft           !   -      -
83      REAL(dp) ::          zavfs    !   -      -
[4990]84      REAL(wp) ::   zavdt, zavds    !   -      -
[9019]85      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
[3]86      !!----------------------------------------------------------------------
[3294]87      !
[3]88      !                                                ! ===============
89      DO jk = 2, jpkm1                                 ! Horizontal slab
90         !                                             ! ===============
91         ! Define the mask
92         ! ---------------
[9019]93!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
94!!gm                     ==>>>  test in the loop instead of use of mask arrays
95!!gm                            and many acces in memory
96         
[13998]97         DO_2D( 1, 1, 1, 1 )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
[12377]98            zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
[13237]99!!gm please, use e3w at Kmm below
[12377]100               &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
101            !
102            zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
103                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
104            zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
105                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
106            !
107            zdt = zaw * ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )
108            zds = zbw * ( ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) 
109            IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
110            zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
111         END_2D
[3]112
[13998]113         DO_2D( 1, 1, 1, 1 )           !==  indicators  ==!
[12377]114            ! stability indicator: msks=1 if rn2>0; 0 elsewhere
115            IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
116            ELSE                                       ;   zmsks(ji,jj) = 1._wp
117            ENDIF
118            ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
119            IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
120            ELSE                                       ;   zmskf(ji,jj) = 1._wp
121            ENDIF
122            ! diffusive layering indicators:
123            !     ! mskdl1=1 if 0< R <1; 0 elsewhere
124            IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
125            ELSE                                       ;   zmskd1(ji,jj) = 1._wp
126            ENDIF
127            !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
128            IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
129            ELSE                                       ;   zmskd2(ji,jj) = 1._wp
130            ENDIF
131            !   mskdl3=1 if 0.5< R <1; 0 elsewhere
132            IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
133            ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
134            ENDIF
135         END_2D
[3]136         ! mask zmsk in order to have avt and avs masked
[7753]137         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
[3]138
139
140         ! Update avt and avs
141         ! ------------------
142         ! Constant eddy coefficient: reset to the background value
[13295]143         DO_2D( 1, 1, 1, 1 )
[12377]144            zinr = 1._wp / zrau(ji,jj)
145            ! salt fingering
146            zrr = zrau(ji,jj) / rn_hsbfr
147            zrr = zrr * zrr
148            zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
149            zavft = 0.7 * zavfs * zinr
150            ! diffusive layering
151            zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
152            zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
153               &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
154            ! add to the eddy viscosity coef. previously computed
155            p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
156            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
157            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
158         END_2D
[3]159         !                                                ! ===============
160      END DO                                              !   End of slab
161      !                                                   ! ===============
[1601]162      !
[12377]163      IF(sn_cfctl%l_prtctl) THEN
[9440]164         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
[49]165      ENDIF
[1601]166      !
[3]167   END SUBROUTINE zdf_ddm
168   
169   !!======================================================================
170END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.