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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 2027

Last change on this file since 2027 was 2027, checked in by cetlod, 14 years ago

Reorganisation of the initialisation phase, see ticket:695
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.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   !!----------------------------------------------------------------------
9#if defined key_zdfddm   ||   defined key_esopa
10   !!----------------------------------------------------------------------
11   !!   'key_zdfddm' :                                     double diffusion
12   !!----------------------------------------------------------------------
13   !!   zdf_ddm       : compute the Ks for salinity
14   !!   zdf_ddm_init  : read namelist and control the parameters
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 in_out_manager  ! I/O manager
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21   USE prtctl          ! Print control
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   zdf_ddm       ! called by step.F90
27   PUBLIC   zdf_ddm_init  ! called by opa.F90
28
29   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
30
31   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avs    !: salinity vertical diffusivity coeff. at w-point
32   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rrau   !: heat/salt buoyancy flux ratio
33
34   !                                  !!* Namelist namzdf_ddm : double diffusive mixing *
35   REAL(wp) ::   rn_avts  = 1.e-4_wp   ! maximum value of avs for salt fingering
36   REAL(wp) ::   rn_hsbfr = 1.6_wp     ! heat/salt buoyancy flux ratio
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
42   !! $Id$
43   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE zdf_ddm( kt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE zdf_ddm  ***
51      !!                   
52      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
53      !!              effect of salt fingering and diffusive convection.
54      !!
55      !! ** Method  :   Diapycnal mixing is increased in case of double
56      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
57      !!      following Merryfield et al. (1999). The rate of double diffusive
58      !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S]
59      !!      which is computed in rn2.F
60      !!         * salt fingering (Schmitt 1981):
61      !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )
62      !!      for Rrau > 1 and rn2 > 0 : zavfs = O
63      !!      otherwise                : zavft = 0.7 zavs / Rrau
64      !!         * diffusive layering (Federov 1988):
65      !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6 
66      !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) )
67      !!      otherwise                   : zavdt = 0
68      !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85)
69      !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau     
70      !!      otherwise                     : zavds = 0
71      !!         * update the eddy diffusivity:
72      !!      avt = avt + zavft + zavdt
73      !!      avs = avs + zavfs + zavds
74      !!      avmu, avmv are required to remain at least above avt and avs.
75      !!     
76      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
77      !!
78      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
79      !!----------------------------------------------------------------------
80      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
81      !!
82      INTEGER  ::   ji, jj , jk     ! dummy loop indices
83      REAL(wp) ::   zinr, zrr       ! temporary scalars
84      REAL(wp) ::   zavft, zavfs    !    -         -
85      REAL(wp) ::   zavdt, zavds    !    -         -
86      REAL(wp), DIMENSION(jpi,jpj) ::   zmsks, zmskf, zmskd1, zmskd2, zmskd3   ! 2D workspace
87      !!----------------------------------------------------------------------
88
89      !                                                ! ===============
90      DO jk = 2, jpkm1                                 ! Horizontal slab
91         !                                             ! ===============
92         ! Define the mask
93         ! ---------------
94         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
95
96         DO jj = 1, jpj                                     ! indicators:
97            DO ji = 1, jpi
98               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
99               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN
100                  zmsks(ji,jj) = 0.e0
101               ELSE
102                  zmsks(ji,jj) = 1.e0
103               ENDIF
104               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
105               IF( rrau(ji,jj,jk) <= 1. ) THEN
106                  zmskf(ji,jj) = 0.e0
107               ELSE
108                  zmskf(ji,jj) = 1.e0
109               ENDIF
110               ! diffusive layering indicators:
111               !   mskdl1=1 if 0<rrau<1; 0 elsewhere
112               IF( rrau(ji,jj,jk) >= 1. ) THEN
113                  zmskd1(ji,jj) = 0.e0
114               ELSE
115                  zmskd1(ji,jj) = 1.e0
116               ENDIF
117               !   mskdl2=1 if 0<rrau<0.5; 0 elsewhere
118               IF( rrau(ji,jj,jk) >= 0.5 ) THEN
119                  zmskd2(ji,jj) = 0.e0
120               ELSE
121                  zmskd2(ji,jj) = 1.e0
122               ENDIF
123               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
124               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN
125                  zmskd3(ji,jj) = 0.e0
126               ELSE
127                  zmskd3(ji,jj) = 1.e0
128               ENDIF
129            END DO
130         END DO
131         ! mask zmsk in order to have avt and avs masked
132         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
133
134
135         ! Update avt and avs
136         ! ------------------
137         ! Constant eddy coefficient: reset to the background value
138!CDIR NOVERRCHK
139         DO jj = 1, jpj
140!CDIR NOVERRCHK
141            DO ji = 1, jpi
142               zinr = 1./rrau(ji,jj,jk)
143               ! salt fingering
144               zrr = rrau(ji,jj,jk)/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 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
151                  &                            +  0.15 * rrau(ji,jj,jk)          * 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) )  * umask(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) )  * vmask(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. )        ! Lateral boundary conditions   (unchanged sign)
178      CALL lbc_lnk( avs , 'W', 1. )
179      CALL lbc_lnk( avm , 'W', 1. )
180      CALL lbc_lnk( avmu, 'U', 1. ) 
181      CALL lbc_lnk( avmv, 'V', 1. )
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   END SUBROUTINE zdf_ddm
190   
191   
192   SUBROUTINE zdf_ddm_init
193      !!----------------------------------------------------------------------
194      !!                  ***  ROUTINE zdf_ddm_init  ***
195      !!
196      !! ** Purpose :   Initialization of double diffusion mixing scheme
197      !!
198      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
199      !!              called by zdf_ddm at the first timestep (nit000)
200      !!----------------------------------------------------------------------
201      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
202      !!----------------------------------------------------------------------
203      !
204      REWIND ( numnam )               ! Read Namelist namzdf_ddm : double diffusion mixing scheme
205      READ   ( numnam, namzdf_ddm )
206      !
207      IF(lwp) THEN                    ! Parameter print
208         WRITE(numout,*)
209         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
210         WRITE(numout,*) '~~~~~~~'
211         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
212         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
213         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
214         WRITE(numout,*)
215      ENDIF
216      !
217   END SUBROUTINE zdf_ddm_init
218
219#else
220   !!----------------------------------------------------------------------
221   !!   Default option :          Dummy module          No double diffusion
222   !!----------------------------------------------------------------------
223   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
224CONTAINS
225   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
226      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
227   END SUBROUTINE zdf_ddm
228#endif
229
230   !!======================================================================
231END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.