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 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

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