source: NEMO/trunk/src/OCE/ZDF/zdfddm.F90 @ 13226

Last change on this file since 13226 was 13226, checked in by orioltp, 3 months ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk

  • Property svn:keywords set to Id
File size: 8.7 KB
Line 
1MODULE zdfddm
2   !!======================================================================
3   !!                       ***  MODULE  zdfddm  ***
4   !! Ocean physics : double diffusion mixing parameterization
5   !!======================================================================
6   !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
8   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
9   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta
10   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   zdf_ddm       : compute the Kz for salinity
15   !!----------------------------------------------------------------------
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
19   USE eosbn2         ! equation of state
20   !
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
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   zdf_ddm       ! called by step.F90
30
31   !! * Substitutions
32#  include "do_loop_substitute.h90"
33   !!----------------------------------------------------------------------
34   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
35   !! $Id$
36   !! Software governed by the CeCILL license (see ./LICENSE)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE zdf_ddm( kt, Kmm, p_avm, p_avt, p_avs )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE zdf_ddm  ***
43      !!                   
44      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
45      !!              effect of salt fingering and diffusive convection.
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
50      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
51      !!         * salt fingering (Schmitt 1981):
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
55      !!         * diffusive layering (Federov 1988):
56      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
57      !!      otherwise                   : zavdt = 0
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     
60      !!      otherwise                     : zavds = 0
61      !!         * update the eddy diffusivity:
62      !!      avt = avt + zavft + zavdt
63      !!      avs = avs + zavfs + zavds
64      !!      avm is required to remain at least above avt and avs.
65      !!     
66      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
67      !!
68      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
69      !!----------------------------------------------------------------------
70      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
71      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
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)
75      !
76      INTEGER  ::   ji, jj , jk     ! dummy loop indices
77      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
78      REAL(wp) ::   zdt, zds
79      REAL(wp) ::   zinr            !   -      -
80      REAL(dp) ::         zrr       !   -      -
81      REAL(wp) ::   zavft           !   -      -
82      REAL(dp) ::          zavfs    !   -      -
83      REAL(wp) ::   zavdt, zavds    !   -      -
84      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
85      !!----------------------------------------------------------------------
86      !
87      !                                                ! ===============
88      DO jk = 2, jpkm1                                 ! Horizontal slab
89         !                                             ! ===============
90         ! Define the mask
91         ! ---------------
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         
96         DO_2D_11_11
97            zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
98!!gm please, use e3w(:,:,:,Kmm) below
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
111
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
135         ! mask zmsk in order to have avt and avs masked
136         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
137
138
139         ! Update avt and avs
140         ! ------------------
141         ! Constant eddy coefficient: reset to the background value
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
158         !                                                ! ===============
159      END DO                                              !   End of slab
160      !                                                   ! ===============
161      !
162      IF(sn_cfctl%l_prtctl) THEN
163         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
164      ENDIF
165      !
166   END SUBROUTINE zdf_ddm
167   
168   !!======================================================================
169END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.