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/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

#1260 : see the associated wiki page for explanation
  • Property svn:keywords set to Id
File size: 12.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   !!----------------------------------------------------------------------
10#if defined key_zdfddm   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_zdfddm' :                                     double diffusion
13   !!----------------------------------------------------------------------
14   !!   zdf_ddm       : compute the Ks for salinity
15   !!   zdf_ddm_init  : read namelist and control the parameters
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers variables
18   USE dom_oce         ! ocean space and time domain variables
19   USE zdf_oce         ! ocean vertical physics variables
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   PUBLIC   zdf_ddm_init  ! called by opa.F90
32   PUBLIC   zdf_ddm_alloc ! called by nemogcm.F90
33
34   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
35
36   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs    !: salinity vertical diffusivity coeff. at w-point
37   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   rrau   !: heat/salt buoyancy flux ratio
38
39   !                      !!* Namelist namzdf_ddm : double diffusive mixing *
40   REAL(wp) ::   rn_avts   ! maximum value of avs for salt fingering
41   REAL(wp) ::   rn_hsbfr  ! heat/salt buoyancy flux ratio
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   INTEGER FUNCTION zdf_ddm_alloc()
53      !!----------------------------------------------------------------------
54      !!                ***  ROUTINE zdf_ddm_alloc  ***
55      !!----------------------------------------------------------------------
56      ALLOCATE( avs(jpi,jpj,jpk), rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc )
57      !
58      IF( lk_mpp             )   CALL mpp_sum ( zdf_ddm_alloc )
59      IF( zdf_ddm_alloc /= 0 )   CALL ctl_warn('zdf_ddm_alloc: failed to allocate arrays')
60   END FUNCTION zdf_ddm_alloc
61
62
63   SUBROUTINE zdf_ddm( kt )
64      !!----------------------------------------------------------------------
65      !!                  ***  ROUTINE zdf_ddm  ***
66      !!                   
67      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
68      !!              effect of salt fingering and diffusive convection.
69      !!
70      !! ** Method  :   Diapycnal mixing is increased in case of double
71      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
72      !!      following Merryfield et al. (1999). The rate of double diffusive
73      !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S]
74      !!      which is computed in rn2.F
75      !!         * salt fingering (Schmitt 1981):
76      !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )
77      !!      for Rrau > 1 and rn2 > 0 : zavfs = O
78      !!      otherwise                : zavft = 0.7 zavs / Rrau
79      !!         * diffusive layering (Federov 1988):
80      !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6 
81      !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) )
82      !!      otherwise                   : zavdt = 0
83      !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85)
84      !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau     
85      !!      otherwise                     : zavds = 0
86      !!         * update the eddy diffusivity:
87      !!      avt = avt + zavft + zavdt
88      !!      avs = avs + zavfs + zavds
89      !!      avmu, avmv are required to remain at least above avt and avs.
90      !!     
91      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
92      !!
93      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
94      !!----------------------------------------------------------------------
95      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
96      !
97      INTEGER  ::   ji, jj , jk     ! dummy loop indices
98      REAL(wp) ::   zinr, zrr       ! temporary scalars
99      REAL(wp) ::   zavft, zavfs    !    -         -
100      REAL(wp) ::   zavdt, zavds    !    -         -
101      REAL(wp), POINTER, DIMENSION(:,:) ::   zmsks, zmskf, zmskd1, zmskd2, zmskd3
102      !!----------------------------------------------------------------------
103      !
104      IF( nn_timing == 1 )  CALL timing_start('zdf_ddm')
105      !
106      CALL wrk_alloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
107
108      !                                                ! ===============
109      DO jk = 2, jpkm1                                 ! Horizontal slab
110         !                                             ! ===============
111         ! Define the mask
112         ! ---------------
113         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
114
115         DO jj = 1, jpj                                     ! indicators:
116            DO ji = 1, jpi
117               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
118               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
119               ELSE                                       ;   zmsks(ji,jj) = 1._wp
120               ENDIF
121               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
122               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp
123               ELSE                                       ;   zmskf(ji,jj) = 1._wp
124               ENDIF
125               ! diffusive layering indicators:
126               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
127               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp
128               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
129               ENDIF
130               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
131               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp
132               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
133               ENDIF
134               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
135               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
136               ELSE                                                         ;   zmskd3(ji,jj) = 1._wp
137               ENDIF
138            END DO
139         END DO
140         ! mask zmsk in order to have avt and avs masked
141         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
142
143
144         ! Update avt and avs
145         ! ------------------
146         ! Constant eddy coefficient: reset to the background value
147         DO jj = 1, jpj
148            DO ji = 1, jpi
149               zinr = 1./rrau(ji,jj,jk)
150               ! salt fingering
151               zrr = rrau(ji,jj,jk)/rn_hsbfr
152               zrr = zrr * zrr
153               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
154               zavft = 0.7 * zavfs * zinr
155               ! diffusive layering
156               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
157               zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
158                  &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
159               ! add to the eddy viscosity coef. previously computed
160               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
161               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
162               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
163            END DO
164         END DO
165
166
167         ! Increase avmu, avmv if necessary
168         ! --------------------------------
169!!gm to be changed following the definition of avm.
170         DO jj = 1, jpjm1
171            DO ji = 1, fs_jpim1   ! vector opt.
172               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
173                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
174                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
175               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
176                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
177                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk)
178            END DO
179         END DO
180         !                                                ! ===============
181      END DO                                              !   End of slab
182      !                                                   ! ===============
183      !
184      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
185      CALL lbc_lnk( avs , 'W', 1._wp )
186      CALL lbc_lnk( avm , 'W', 1._wp )
187      CALL lbc_lnk( avmu, 'U', 1._wp ) 
188      CALL lbc_lnk( avmv, 'V', 1._wp )
189
190      IF(ln_ctl) THEN
191         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
192         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
193            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
194      ENDIF
195      !
196      CALL wrk_dealloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
197      !
198      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
199      !
200   END SUBROUTINE zdf_ddm
201   
202   
203   SUBROUTINE zdf_ddm_init
204      !!----------------------------------------------------------------------
205      !!                  ***  ROUTINE zdf_ddm_init  ***
206      !!
207      !! ** Purpose :   Initialization of double diffusion mixing scheme
208      !!
209      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
210      !!              called by zdf_ddm at the first timestep (nit000)
211      !!----------------------------------------------------------------------
212      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
213      INTEGER  ::   ios                 ! Local integer output status for namelist read
214      !!----------------------------------------------------------------------
215      !
216      REWIND( numnam_ref )              ! Namelist namzdf_ddm in reference namelist : Double diffusion mixing scheme
217      READ  ( numnam_ref, namzdf_ddm, IOSTAT = ios, ERR = 901)
218901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ddm in reference namelist', lwp )
219
220      REWIND( numnam_cfg )              ! Namelist namzdf_ddm in configuration namelist : Double diffusion mixing scheme
221      READ  ( numnam_cfg, namzdf_ddm, IOSTAT = ios, ERR = 902 )
222902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ddm in configuration namelist', lwp )
223      WRITE ( numond, namzdf_ddm )
224      !
225      IF(lwp) THEN                    ! Parameter print
226         WRITE(numout,*)
227         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
228         WRITE(numout,*) '~~~~~~~'
229         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
230         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
231         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
232      ENDIF
233      !
234      !                               ! allocate zdfddm arrays
235      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
236      !                               ! initialization to masked Kz
237      avs(:,:,:) = rn_avt0 * tmask(:,:,:) 
238      !
239   END SUBROUTINE zdf_ddm_init
240
241#else
242   !!----------------------------------------------------------------------
243   !!   Default option :          Dummy module          No double diffusion
244   !!----------------------------------------------------------------------
245   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
246CONTAINS
247   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
248      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
249   END SUBROUTINE zdf_ddm
250   SUBROUTINE zdf_ddm_init            ! Dummy routine
251   END SUBROUTINE zdf_ddm_init
252#endif
253
254   !!======================================================================
255END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.