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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

set proper svn properties to all files...
  • Property svn:keywords set to Id
File size: 11.0 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
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   zdf_ddm       ! called by step.F90
28   PUBLIC   zdf_ddm_init  ! called by opa.F90
29
30   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
31
32   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avs    !: salinity vertical diffusivity coeff. at w-point
33   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rrau   !: heat/salt buoyancy flux ratio
34
35   !                                  !!* Namelist namzdf_ddm : double diffusive mixing *
36   REAL(wp) ::   rn_avts  = 1.e-4_wp   ! maximum value of avs for salt fingering
37   REAL(wp) ::   rn_hsbfr = 1.6_wp     ! heat/salt buoyancy flux ratio
38
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3,3 , LOCEAN-IPSL (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE zdf_ddm( kt )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE zdf_ddm  ***
52      !!                   
53      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
54      !!              effect of salt fingering and diffusive convection.
55      !!
56      !! ** Method  :   Diapycnal mixing is increased in case of double
57      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
58      !!      following Merryfield et al. (1999). The rate of double diffusive
59      !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S]
60      !!      which is computed in rn2.F
61      !!         * salt fingering (Schmitt 1981):
62      !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )
63      !!      for Rrau > 1 and rn2 > 0 : zavfs = O
64      !!      otherwise                : zavft = 0.7 zavs / Rrau
65      !!         * diffusive layering (Federov 1988):
66      !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6 
67      !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) )
68      !!      otherwise                   : zavdt = 0
69      !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85)
70      !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau     
71      !!      otherwise                     : zavds = 0
72      !!         * update the eddy diffusivity:
73      !!      avt = avt + zavft + zavdt
74      !!      avs = avs + zavfs + zavds
75      !!      avmu, avmv are required to remain at least above avt and avs.
76      !!     
77      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
78      !!
79      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
80      !!----------------------------------------------------------------------
81      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
82      !!
83      INTEGER  ::   ji, jj , jk     ! dummy loop indices
84      REAL(wp) ::   zinr, zrr       ! temporary scalars
85      REAL(wp) ::   zavft, zavfs    !    -         -
86      REAL(wp) ::   zavdt, zavds    !    -         -
87      REAL(wp), DIMENSION(jpi,jpj) ::   zmsks, zmskf, zmskd1, zmskd2, zmskd3   ! 2D workspace
88      !!----------------------------------------------------------------------
89
90      !                                                ! ===============
91      DO jk = 2, jpkm1                                 ! Horizontal slab
92         !                                             ! ===============
93         ! Define the mask
94         ! ---------------
95         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
96
97         DO jj = 1, jpj                                     ! indicators:
98            DO ji = 1, jpi
99               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
100               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0.e0
101               ELSE                                       ;   zmsks(ji,jj) = 1.e0
102               ENDIF
103               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
104               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0.e0
105               ELSE                                       ;   zmskf(ji,jj) = 1.e0
106               ENDIF
107               ! diffusive layering indicators:
108               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
109               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0.e0
110               ELSE                                       ;   zmskd1(ji,jj) = 1.e0
111               ENDIF
112               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
113               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0.e0
114               ELSE                                       ;   zmskd2(ji,jj) = 1.e0
115               ENDIF
116               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
117               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0.e0
118               ELSE                                                         ;   zmskd3(ji,jj) = 1.e0
119               ENDIF
120            END DO
121         END DO
122         ! mask zmsk in order to have avt and avs masked
123         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
124
125
126         ! Update avt and avs
127         ! ------------------
128         ! Constant eddy coefficient: reset to the background value
129!CDIR NOVERRCHK
130         DO jj = 1, jpj
131!CDIR NOVERRCHK
132            DO ji = 1, jpi
133               zinr = 1./rrau(ji,jj,jk)
134               ! salt fingering
135               zrr = rrau(ji,jj,jk)/rn_hsbfr
136               zrr = zrr * zrr
137               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
138               zavft = 0.7 * zavfs * zinr
139               ! diffusive layering
140               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
141               zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
142                  &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
143               ! add to the eddy viscosity coef. previously computed
144               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
145               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
146               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
147            END DO
148         END DO
149
150
151         ! Increase avmu, avmv if necessary
152         ! --------------------------------
153!!gm to be changed following the definition of avm.
154         DO jj = 1, jpjm1
155            DO ji = 1, fs_jpim1   ! vector opt.
156               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
157                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
158                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
159               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
160                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
161                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk)
162            END DO
163         END DO
164         !                                                ! ===============
165      END DO                                              !   End of slab
166      !                                                   ! ===============
167      !
168      CALL lbc_lnk( avt , 'W', 1. )        ! Lateral boundary conditions   (unchanged sign)
169      CALL lbc_lnk( avs , 'W', 1. )
170      CALL lbc_lnk( avm , 'W', 1. )
171      CALL lbc_lnk( avmu, 'U', 1. ) 
172      CALL lbc_lnk( avmv, 'V', 1. )
173
174      IF(ln_ctl) THEN
175         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
176         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
177            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
178      ENDIF
179      !
180   END SUBROUTINE zdf_ddm
181   
182   
183   SUBROUTINE zdf_ddm_init
184      !!----------------------------------------------------------------------
185      !!                  ***  ROUTINE zdf_ddm_init  ***
186      !!
187      !! ** Purpose :   Initialization of double diffusion mixing scheme
188      !!
189      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
190      !!              called by zdf_ddm at the first timestep (nit000)
191      !!----------------------------------------------------------------------
192      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
193      !!----------------------------------------------------------------------
194      !
195      REWIND ( numnam )               ! Read Namelist namzdf_ddm : double diffusion mixing scheme
196      READ   ( numnam, namzdf_ddm )
197      !
198      IF(lwp) THEN                    ! Parameter print
199         WRITE(numout,*)
200         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
201         WRITE(numout,*) '~~~~~~~'
202         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
203         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
204         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
205         WRITE(numout,*)
206      ENDIF
207      !
208   END SUBROUTINE zdf_ddm_init
209
210#else
211   !!----------------------------------------------------------------------
212   !!   Default option :          Dummy module          No double diffusion
213   !!----------------------------------------------------------------------
214   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
215CONTAINS
216   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
217      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
218   END SUBROUTINE zdf_ddm
219   SUBROUTINE zdf_ddm_init            ! Dummy routine
220   END SUBROUTINE zdf_ddm_init
221#endif
222
223   !!======================================================================
224END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.