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 branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 7931

Last change on this file since 7931 was 7931, checked in by gm, 7 years ago

#1880 (HPC-09): remove key_zdfddm + phasing with last changes of HPC08 branch
  • Property svn:keywords set to Id
File size: 9.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   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   zdf_ddm       : compute the Kz for salinity
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers variables
16   USE dom_oce        ! ocean space and time domain variables
17   USE zdf_oce        ! ocean vertical physics variables
18   USE eosbn2         ! equation of state
19   !
20   USE in_out_manager ! I/O manager
21   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
22   USE prtctl         ! Print control
23   USE lib_mpp        ! MPP library
24   USE wrk_nemo       ! work arrays
25   USE timing         ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   zdf_ddm       ! called by step.F90
31
32   !! * Substitutions
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
36   !! $Id$
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE zdf_ddm( kt )
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE zdf_ddm  ***
44      !!                   
45      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
46      !!              effect of salt fingering and diffusive convection.
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
51      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
52      !!         * salt fingering (Schmitt 1981):
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
56      !!         * diffusive layering (Federov 1988):
57      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
58      !!      otherwise                   : zavdt = 0
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     
61      !!      otherwise                     : zavds = 0
62      !!         * update the eddy diffusivity:
63      !!      avt = avt + zavft + zavdt
64      !!      avs = avs + zavfs + zavds
65      !!      avmu, avmv are required to remain at least above avt and avs.
66      !!     
67      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
68      !!
69      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
70      !!----------------------------------------------------------------------
71      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
72      !
73      INTEGER  ::   ji, jj , jk     ! dummy loop indices
74      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
75      REAL(wp) ::   zdt, zds
76      REAL(wp) ::   zinr, zrr       !   -      -
77      REAL(wp) ::   zavft, zavfs    !   -      -
78      REAL(wp) ::   zavdt, zavds    !   -      -
79      REAL(wp), POINTER, DIMENSION(:,:) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
80      !!----------------------------------------------------------------------
81      !
82      IF( nn_timing == 1 )  CALL timing_start('zdf_ddm')
83      !
84      CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
85      !
86      !                                                ! ===============
87      DO jk = 2, jpkm1                                 ! Horizontal slab
88         !                                             ! ===============
89         ! Define the mask
90         ! ---------------
91         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s])
92            DO ji = 1, jpi
93               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
94                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
95               !
96               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
97                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
98               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
99                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
100               !
101               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
102               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
103               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
104               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
105            END DO
106         END DO
107
108         DO jj = 1, jpj                                     ! indicators:
109            DO ji = 1, jpi
110               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
111               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
112               ELSE                                       ;   zmsks(ji,jj) = 1._wp
113               ENDIF
114               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
115               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
116               ELSE                                       ;   zmskf(ji,jj) = 1._wp
117               ENDIF
118               ! diffusive layering indicators:
119               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
120               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
121               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
122               ENDIF
123               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
124               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
125               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
126               ENDIF
127               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
128               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
129               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
130               ENDIF
131            END DO
132         END DO
133         ! mask zmsk in order to have avt and avs masked
134         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
135
136
137         ! Update avt and avs
138         ! ------------------
139         ! Constant eddy coefficient: reset to the background value
140         DO jj = 1, jpj
141            DO ji = 1, jpi
142               zinr = 1._wp / zrau(ji,jj)
143               ! salt fingering
144               zrr = zrau(ji,jj) / rn_hsbfr
145               zrr = zrr * zrr
146               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
147               zavft = 0.7 * zavfs * zinr
148               ! diffusive layering
149               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
150               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
151                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
152               ! add to the eddy viscosity coef. previously computed
153               avs(ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
154               avt(ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
155               avm(ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
156            END DO
157         END DO
158
159
160         ! Increase avmu, avmv if necessary
161         ! --------------------------------
162!!gm to be changed following the definition of avm.
163         DO jj = 1, jpjm1
164            DO ji = 1, fs_jpim1   ! vector opt.
165               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
166                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
167                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * wumask(ji,jj,jk)
168               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
169                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
170                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * wvmask(ji,jj,jk)
171            END DO
172         END DO
173         !                                                ! ===============
174      END DO                                              !   End of slab
175      !                                                   ! ===============
176      !
177      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
178      CALL lbc_lnk( avs , 'W', 1._wp )
179      CALL lbc_lnk( avm , 'W', 1._wp )
180      CALL lbc_lnk( avmu, 'U', 1._wp ) 
181      CALL lbc_lnk( avmv, 'V', 1._wp )
182
183      IF(ln_ctl) THEN
184         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
185         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
186            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
187      ENDIF
188      !
189      CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
190      !
191      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
192      !
193   END SUBROUTINE zdf_ddm
194   
195   !!======================================================================
196END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.