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

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 4704

Last change on this file since 4704 was 4666, checked in by mathiot, 10 years ago

#1331 : add ISOMIP config files + ice shelf code
  • Property svn:keywords set to Id
File size: 12.9 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 4.0 , NEMO Consortium (2011)
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      !                                                ! Horizontal slab
110         !                                             ! ===============
111      DO jj = 1, jpj                                     ! indicators:
112         DO ji = 1, jpi
113            DO jk = mikt(ji,jj)+1, jpkm1                                 ! Horizontal slab
114         ! Define the mask
115         ! ---------------
116               rrau(ji,jj,jk) = MAX( 1.e-20, rrau(ji,jj,jk) )         ! only retains positive value of rrau
117
118               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
119               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
120               ELSE                                       ;   zmsks(ji,jj) = 1._wp
121               ENDIF
122               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
123               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp
124               ELSE                                       ;   zmskf(ji,jj) = 1._wp
125               ENDIF
126               ! diffusive layering indicators:
127               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
128               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp
129               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
130               ENDIF
131               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
132               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp
133               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
134               ENDIF
135               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
136               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
137               ELSE                                                         ;   zmskd3(ji,jj) = 1._wp
138               ENDIF
139         ! mask zmsk in order to have avt and avs masked
140               zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
141
142
143         ! Update avt and avs
144         ! ------------------
145         ! Constant eddy coefficient: reset to the background value
146               zinr = 1./rrau(ji,jj,jk)
147               ! salt fingering
148               zrr = rrau(ji,jj,jk)/rn_hsbfr
149               zrr = zrr * zrr
150               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
151               zavft = 0.7 * zavfs * zinr
152               ! diffusive layering
153               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
154               zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
155                  &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
156               ! add to the eddy viscosity coef. previously computed
157               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
158               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
159               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
160            END DO
161         END DO
162      END DO
163         ! Increase avmu, avmv if necessary
164         ! --------------------------------
165!!gm to be changed following the definition of avm.
166      DO jj = 1, jpjm1
167         DO ji = 1, fs_jpim1   ! vector opt.
168            DO jk = miku(ji,jj)+1, jpkm1                                 ! Horizontal slab
169               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
170                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
171                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
172            END DO
173            DO jk = mikv(ji,jj)+1, jpkm1
174               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
175                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
176                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk)
177            END DO
178         END DO
179         !                                                ! ===============
180      END DO                                              !   End of slab
181      !                                                   ! ===============
182      !
183      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
184      CALL lbc_lnk( avs , 'W', 1._wp )
185      CALL lbc_lnk( avm , 'W', 1._wp )
186      CALL lbc_lnk( avmu, 'U', 1._wp ) 
187      CALL lbc_lnk( avmv, 'V', 1._wp )
188
189      IF(ln_ctl) THEN
190         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
191         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
192            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
193      ENDIF
194      !
195      CALL wrk_dealloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
196      !
197      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
198      !
199   END SUBROUTINE zdf_ddm
200   
201   
202   SUBROUTINE zdf_ddm_init
203      !!----------------------------------------------------------------------
204      !!                  ***  ROUTINE zdf_ddm_init  ***
205      !!
206      !! ** Purpose :   Initialization of double diffusion mixing scheme
207      !!
208      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
209      !!              called by zdf_ddm at the first timestep (nit000)
210      !!----------------------------------------------------------------------
211      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
212      INTEGER  ::   ios                 ! Local integer output status for namelist read
213      !!----------------------------------------------------------------------
214      !
215      REWIND( numnam_ref )              ! Namelist namzdf_ddm in reference namelist : Double diffusion mixing scheme
216      READ  ( numnam_ref, namzdf_ddm, IOSTAT = ios, ERR = 901)
217901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ddm in reference namelist', lwp )
218
219      REWIND( numnam_cfg )              ! Namelist namzdf_ddm in configuration namelist : Double diffusion mixing scheme
220      READ  ( numnam_cfg, namzdf_ddm, IOSTAT = ios, ERR = 902 )
221902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ddm in configuration namelist', lwp )
222      IF(lwm) WRITE ( numond, namzdf_ddm )
223      !
224      IF(lwp) THEN                    ! Parameter print
225         WRITE(numout,*)
226         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
227         WRITE(numout,*) '~~~~~~~'
228         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
229         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
230         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
231      ENDIF
232      !
233      !                               ! allocate zdfddm arrays
234      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
235      !                               ! initialization to masked Kz
236      avs(:,:,:) = rn_avt0 * tmask(:,:,:) 
237      !
238   END SUBROUTINE zdf_ddm_init
239
240#else
241   !!----------------------------------------------------------------------
242   !!   Default option :          Dummy module          No double diffusion
243   !!----------------------------------------------------------------------
244   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
245CONTAINS
246   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
247      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
248   END SUBROUTINE zdf_ddm
249   SUBROUTINE zdf_ddm_init            ! Dummy routine
250   END SUBROUTINE zdf_ddm_init
251#endif
252
253   !!======================================================================
254END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.