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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

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