source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 3229

Last change on this file since 3229 was 3229, checked in by charris, 10 years ago

Added timing calls to most significant routines in LDF, SBC and ZDF.

  • Property svn:keywords set to Id
File size: 12.1 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  = 1.e-4_wp   ! maximum value of avs for salt fingering
41   REAL(wp) ::   rn_hsbfr = 1.6_wp     ! 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      DO jk = 2, jpkm1                                 ! Horizontal slab
110         !                                             ! ===============
111         ! Define the mask
112         ! ---------------
113         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
114
115         DO jj = 1, jpj                                     ! indicators:
116            DO ji = 1, jpi
117               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
118               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
119               ELSE                                       ;   zmsks(ji,jj) = 1._wp
120               ENDIF
121               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
122               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp
123               ELSE                                       ;   zmskf(ji,jj) = 1._wp
124               ENDIF
125               ! diffusive layering indicators:
126               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
127               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp
128               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
129               ENDIF
130               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
131               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp
132               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
133               ENDIF
134               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
135               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
136               ELSE                                                         ;   zmskd3(ji,jj) = 1._wp
137               ENDIF
138            END DO
139         END DO
140         ! mask zmsk in order to have avt and avs masked
141         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
142
143
144         ! Update avt and avs
145         ! ------------------
146         ! Constant eddy coefficient: reset to the background value
147!CDIR NOVERRCHK
148         DO jj = 1, jpj
149!CDIR NOVERRCHK
150            DO ji = 1, jpi
151               zinr = 1./rrau(ji,jj,jk)
152               ! salt fingering
153               zrr = rrau(ji,jj,jk)/rn_hsbfr
154               zrr = zrr * zrr
155               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
156               zavft = 0.7 * zavfs * zinr
157               ! diffusive layering
158               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
159               zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
160                  &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
161               ! add to the eddy viscosity coef. previously computed
162               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
163               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
164               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
165            END DO
166         END DO
167
168
169         ! Increase avmu, avmv if necessary
170         ! --------------------------------
171!!gm to be changed following the definition of avm.
172         DO jj = 1, jpjm1
173            DO ji = 1, fs_jpim1   ! vector opt.
174               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
175                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
176                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
177               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
178                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
179                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk)
180            END DO
181         END DO
182         !                                                ! ===============
183      END DO                                              !   End of slab
184      !                                                   ! ===============
185      !
186      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
187      CALL lbc_lnk( avs , 'W', 1._wp )
188      CALL lbc_lnk( avm , 'W', 1._wp )
189      CALL lbc_lnk( avmu, 'U', 1._wp ) 
190      CALL lbc_lnk( avmv, 'V', 1._wp )
191
192      IF(ln_ctl) THEN
193         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
194         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
195            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
196      ENDIF
197      !
198      CALL wrk_dealloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
199      !
200      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
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      ENDIF
228      !
229      !                              ! allocate zdfddm arrays
230      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
231      !
232   END SUBROUTINE zdf_ddm_init
233
234#else
235   !!----------------------------------------------------------------------
236   !!   Default option :          Dummy module          No double diffusion
237   !!----------------------------------------------------------------------
238   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
239CONTAINS
240   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
241      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
242   END SUBROUTINE zdf_ddm
243   SUBROUTINE zdf_ddm_init            ! Dummy routine
244   END SUBROUTINE zdf_ddm_init
245#endif
246
247   !!======================================================================
248END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.