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

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 9294

Last change on this file since 9294 was 7753, checked in by mocavero, 7 years ago

Reverting trunk to remove OpenMP

  • Property svn:keywords set to Id
File size: 13.7 KB
RevLine 
[3]1MODULE zdfddm
2   !!======================================================================
3   !!                       ***  MODULE  zdfddm  ***
4   !! Ocean physics : double diffusion mixing parameterization
5   !!======================================================================
[1601]6   !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
[4990]8   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
9   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta
[1601]10   !!----------------------------------------------------------------------
[5836]11#if defined key_zdfddm
[3]12   !!----------------------------------------------------------------------
13   !!   'key_zdfddm' :                                     double diffusion
14   !!----------------------------------------------------------------------
15   !!   zdf_ddm       : compute the Ks for salinity
16   !!   zdf_ddm_init  : read namelist and control the parameters
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE zdf_oce         ! ocean vertical physics variables
[4990]21   USE eosbn2         ! equation of state
22   !
[3]23   USE in_out_manager  ! I/O manager
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[258]25   USE prtctl          ! Print control
[3294]26   USE lib_mpp         ! MPP library
27   USE wrk_nemo        ! work arrays
28   USE timing          ! Timing
[3]29
30   IMPLICIT NONE
31   PRIVATE
32
[2528]33   PUBLIC   zdf_ddm       ! called by step.F90
34   PUBLIC   zdf_ddm_init  ! called by opa.F90
[2715]35   PUBLIC   zdf_ddm_alloc ! called by nemogcm.F90
[3]36
[1537]37   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
[3]38
[4990]39   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs   !: salinity vertical diffusivity coeff. at w-point
[3]40
[4990]41   !                       !!* Namelist namzdf_ddm : double diffusive mixing *
42   REAL(wp) ::   rn_avts    ! maximum value of avs for salt fingering
43   REAL(wp) ::   rn_hsbfr   ! heat/salt buoyancy flux ratio
[1537]44
[3]45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
[4990]48   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
[2528]49   !! $Id$
[2715]50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]51   !!----------------------------------------------------------------------
52CONTAINS
53
[2715]54   INTEGER FUNCTION zdf_ddm_alloc()
55      !!----------------------------------------------------------------------
56      !!                ***  ROUTINE zdf_ddm_alloc  ***
57      !!----------------------------------------------------------------------
[4990]58      ALLOCATE( avs(jpi,jpj,jpk) , STAT= zdf_ddm_alloc )
[2715]59      IF( lk_mpp             )   CALL mpp_sum ( zdf_ddm_alloc )
60      IF( zdf_ddm_alloc /= 0 )   CALL ctl_warn('zdf_ddm_alloc: failed to allocate arrays')
61   END FUNCTION zdf_ddm_alloc
62
63
[3]64   SUBROUTINE zdf_ddm( kt )
65      !!----------------------------------------------------------------------
66      !!                  ***  ROUTINE zdf_ddm  ***
67      !!                   
68      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
[1601]69      !!              effect of salt fingering and diffusive convection.
[3]70      !!
71      !! ** Method  :   Diapycnal mixing is increased in case of double
72      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
73      !!      following Merryfield et al. (1999). The rate of double diffusive
[4990]74      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
[3]75      !!         * salt fingering (Schmitt 1981):
[4990]76      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
77      !!      for R > 1 and rn2 > 0 : zavfs = O
78      !!      otherwise                : zavft = 0.7 zavs / R
[3]79      !!         * diffusive layering (Federov 1988):
[4990]80      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
[3]81      !!      otherwise                   : zavdt = 0
[4990]82      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
83      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
[3]84      !!      otherwise                     : zavds = 0
85      !!         * update the eddy diffusivity:
86      !!      avt = avt + zavft + zavdt
87      !!      avs = avs + zavfs + zavds
88      !!      avmu, avmv are required to remain at least above avt and avs.
89      !!     
[1601]90      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
[3]91      !!
[1601]92      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
[3]93      !!----------------------------------------------------------------------
[1601]94      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
[2715]95      !
[1601]96      INTEGER  ::   ji, jj , jk     ! dummy loop indices
[4990]97      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
98      REAL(wp) ::   zdt, zds
99      REAL(wp) ::   zinr, zrr       !   -      -
100      REAL(wp) ::   zavft, zavfs    !   -      -
101      REAL(wp) ::   zavdt, zavds    !   -      -
102      REAL(wp), POINTER, DIMENSION(:,:) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
[3]103      !!----------------------------------------------------------------------
[3294]104      !
105      IF( nn_timing == 1 )  CALL timing_start('zdf_ddm')
106      !
[4990]107      CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
108      !
[3]109      !                                                ! ===============
110      DO jk = 2, jpkm1                                 ! Horizontal slab
111         !                                             ! ===============
112         ! Define the mask
113         ! ---------------
[4990]114         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s])
115            DO ji = 1, jpi
[6140]116               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
117                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
[4990]118               !
119               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
120                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
121               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
122                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
123               !
124               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
125               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
126               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
127               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
128            END DO
129         END DO
[3]130
[1601]131         DO jj = 1, jpj                                     ! indicators:
[3]132            DO ji = 1, jpi
133               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
[2715]134               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
135               ELSE                                       ;   zmsks(ji,jj) = 1._wp
[3]136               ENDIF
[4990]137               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
138               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
[2715]139               ELSE                                       ;   zmskf(ji,jj) = 1._wp
[3]140               ENDIF
141               ! diffusive layering indicators:
[4990]142               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
143               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
[2715]144               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
[3]145               ENDIF
[4990]146               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
147               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
[2715]148               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
[3]149               ENDIF
[4990]150               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
151               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
152               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
[3]153               ENDIF
154            END DO
155         END DO
156         ! mask zmsk in order to have avt and avs masked
[7753]157         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
[3]158
159
160         ! Update avt and avs
161         ! ------------------
162         ! Constant eddy coefficient: reset to the background value
163         DO jj = 1, jpj
164            DO ji = 1, jpi
[4990]165               zinr = 1._wp / zrau(ji,jj)
[3]166               ! salt fingering
[4990]167               zrr = zrau(ji,jj) / rn_hsbfr
[3]168               zrr = zrr * zrr
[1601]169               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
[1163]170               zavft = 0.7 * zavfs * zinr
[3]171               ! diffusive layering
[1601]172               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
[4990]173               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
174                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
[3]175               ! add to the eddy viscosity coef. previously computed
[6497]176# if defined key_zdftmx_new
177               ! key_zdftmx_new: New internal wave-driven param: use avs value computed by zdftmx
178               avs (ji,jj,jk) = avs(ji,jj,jk) + zavfs + zavds
179# else
[3]180               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
[6497]181# endif
[3]182               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
[1527]183               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
[3]184            END DO
185         END DO
186
187
188         ! Increase avmu, avmv if necessary
189         ! --------------------------------
[1527]190!!gm to be changed following the definition of avm.
[3]191         DO jj = 1, jpjm1
192            DO ji = 1, fs_jpim1   ! vector opt.
193               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
[1601]194                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
[5120]195                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * wumask(ji,jj,jk)
[3]196               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
[1601]197                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
[5120]198                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * wvmask(ji,jj,jk)
[3]199            END DO
200         END DO
201         !                                                ! ===============
202      END DO                                              !   End of slab
203      !                                                   ! ===============
[1601]204      !
[2715]205      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
206      CALL lbc_lnk( avs , 'W', 1._wp )
207      CALL lbc_lnk( avm , 'W', 1._wp )
208      CALL lbc_lnk( avmu, 'U', 1._wp ) 
209      CALL lbc_lnk( avmv, 'V', 1._wp )
[49]210
[258]211      IF(ln_ctl) THEN
212         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
[516]213         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
214            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
[49]215      ENDIF
[1601]216      !
[4990]217      CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
[2715]218      !
[3294]219      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
220      !
[3]221   END SUBROUTINE zdf_ddm
222   
223   
224   SUBROUTINE zdf_ddm_init
225      !!----------------------------------------------------------------------
226      !!                  ***  ROUTINE zdf_ddm_init  ***
227      !!
228      !! ** Purpose :   Initialization of double diffusion mixing scheme
229      !!
[1601]230      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
[1537]231      !!              called by zdf_ddm at the first timestep (nit000)
[3]232      !!----------------------------------------------------------------------
[4990]233      INTEGER ::   ios   ! local integer
234      !!
[1601]235      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
[3]236      !!----------------------------------------------------------------------
[1537]237      !
[4147]238      REWIND( numnam_ref )              ! Namelist namzdf_ddm in reference namelist : Double diffusion mixing scheme
239      READ  ( numnam_ref, namzdf_ddm, IOSTAT = ios, ERR = 901)
240901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ddm in reference namelist', lwp )
241
242      REWIND( numnam_cfg )              ! Namelist namzdf_ddm in configuration namelist : Double diffusion mixing scheme
243      READ  ( numnam_cfg, namzdf_ddm, IOSTAT = ios, ERR = 902 )
244902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ddm in configuration namelist', lwp )
[4624]245      IF(lwm) WRITE ( numond, namzdf_ddm )
[1537]246      !
247      IF(lwp) THEN                    ! Parameter print
[3]248         WRITE(numout,*)
249         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
250         WRITE(numout,*) '~~~~~~~'
[1601]251         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
[1537]252         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
253         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
[3]254      ENDIF
[1537]255      !
[3610]256      !                               ! allocate zdfddm arrays
[2715]257      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
[3610]258      !                               ! initialization to masked Kz
[7753]259      avs(:,:,:) = rn_avt0 * wmask(:,:,:) 
[2715]260      !
[3]261   END SUBROUTINE zdf_ddm_init
262
263#else
264   !!----------------------------------------------------------------------
265   !!   Default option :          Dummy module          No double diffusion
266   !!----------------------------------------------------------------------
[16]267   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
[3]268CONTAINS
269   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
[16]270      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
[3]271   END SUBROUTINE zdf_ddm
[2528]272   SUBROUTINE zdf_ddm_init            ! Dummy routine
273   END SUBROUTINE zdf_ddm_init
[3]274#endif
275
276   !!======================================================================
277END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.