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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 2690

Last change on this file since 2690 was 2690, checked in by gm, 13 years ago

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