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

Last change on this file since 3042 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

  • Property svn:keywords set to Id
File size: 12.2 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
[2528]8   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
[1601]9   !!----------------------------------------------------------------------
[3]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)
[258]22   USE prtctl          ! Print control
[2715]23   USE lib_mpp            ! MPP library
[3]24
25   IMPLICIT NONE
26   PRIVATE
27
[2528]28   PUBLIC   zdf_ddm       ! called by step.F90
29   PUBLIC   zdf_ddm_init  ! called by opa.F90
[2715]30   PUBLIC   zdf_ddm_alloc ! called by nemogcm.F90
[3]31
[1537]32   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
[3]33
[2715]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
[3]36
[1601]37   !                                  !!* Namelist namzdf_ddm : double diffusive mixing *
[1537]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
[3]41   !! * Substitutions
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
[2715]44   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[2528]45   !! $Id$
[2715]46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]47   !!----------------------------------------------------------------------
48CONTAINS
49
[2715]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
[3]61   SUBROUTINE zdf_ddm( kt )
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE zdf_ddm  ***
64      !!                   
65      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
[1601]66      !!              effect of salt fingering and diffusive convection.
[3]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):
[1537]74      !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )
[3]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      !!     
[1601]89      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
[3]90      !!
[1601]91      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
[3]92      !!----------------------------------------------------------------------
[2715]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      !
[1601]97      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
[2715]98      !
[1601]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    !    -         -
[3]103      !!----------------------------------------------------------------------
104
[2715]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
[3]109      !                                                ! ===============
110      DO jk = 2, jpkm1                                 ! Horizontal slab
111         !                                             ! ===============
112         ! Define the mask
113         ! ---------------
[1601]114         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
[3]115
[1601]116         DO jj = 1, jpj                                     ! indicators:
[3]117            DO ji = 1, jpi
118               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
[2715]119               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
120               ELSE                                       ;   zmsks(ji,jj) = 1._wp
[3]121               ENDIF
122               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
[2715]123               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp
124               ELSE                                       ;   zmskf(ji,jj) = 1._wp
[3]125               ENDIF
126               ! diffusive layering indicators:
[2528]127               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
[2715]128               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp
129               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
[3]130               ENDIF
[2528]131               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
[2715]132               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp
133               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
[3]134               ENDIF
135               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
[2715]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
[3]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
[1537]154               zrr = rrau(ji,jj,jk)/rn_hsbfr
[3]155               zrr = zrr * zrr
[1601]156               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
[1163]157               zavft = 0.7 * zavfs * zinr
[3]158               ! diffusive layering
[1601]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)  )
[3]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
[1527]165               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
[3]166            END DO
167         END DO
168
169
170         ! Increase avmu, avmv if necessary
171         ! --------------------------------
[1527]172!!gm to be changed following the definition of avm.
[3]173         DO jj = 1, jpjm1
174            DO ji = 1, fs_jpim1   ! vector opt.
175               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
[1601]176                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
177                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
[3]178               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
[1601]179                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
180                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk)
[3]181            END DO
182         END DO
183         !                                                ! ===============
184      END DO                                              !   End of slab
185      !                                                   ! ===============
[1601]186      !
[2715]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 )
[49]192
[258]193      IF(ln_ctl) THEN
194         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
[516]195         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
196            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
[49]197      ENDIF
[1601]198      !
[2715]199      IF( wrk_not_released(2, 1,2,3,4,5) )   CALL ctl_stop('zdf_ddm: Release of workspace arrays failed')
200      !
[3]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      !!
[1601]210      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
[1537]211      !!              called by zdf_ddm at the first timestep (nit000)
[3]212      !!----------------------------------------------------------------------
[1601]213      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
[3]214      !!----------------------------------------------------------------------
[1537]215      !
[2715]216      REWIND( numnam )                ! Read Namelist namzdf_ddm : double diffusion mixing scheme
217      READ  ( numnam, namzdf_ddm )
[1537]218      !
219      IF(lwp) THEN                    ! Parameter print
[3]220         WRITE(numout,*)
221         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
222         WRITE(numout,*) '~~~~~~~'
[1601]223         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
[1537]224         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
225         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
[3]226      ENDIF
[1537]227      !
[2715]228      !                              ! allocate zdfddm arrays
229      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
230      !
[3]231   END SUBROUTINE zdf_ddm_init
232
233#else
234   !!----------------------------------------------------------------------
235   !!   Default option :          Dummy module          No double diffusion
236   !!----------------------------------------------------------------------
[16]237   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
[3]238CONTAINS
239   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
[16]240      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
[3]241   END SUBROUTINE zdf_ddm
[2528]242   SUBROUTINE zdf_ddm_init            ! Dummy routine
243   END SUBROUTINE zdf_ddm_init
[3]244#endif
245
246   !!======================================================================
247END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.