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.
zdfconvadj.F90 in branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfconvadj.F90 @ 6837

Last change on this file since 6837 was 6837, checked in by timgraham, 8 years ago

Added in convective adjustment module and associated changes
File size: 12.0 KB
Line 
1MODULE zdfconvadj
2   !!======================================================================
3   !!                       ***  MODULE  zdfconvadj  ***
4   !! Perform a convective adjustment of an unstable water column at overflow
5   !! grid points.
6   !!======================================================================
7   !! History :  3.4  !  2013-08  (T. Graham)  Original code adapted from UKMO HadCM3 code
8   !!
9   !!----------------------------------------------------------------------
10   USE oce             ! ocean dynamics and tracers variables
11   USE dom_oce         ! ocean space and time domain variables
12   USE zdf_oce         ! ocean vertical physics variables
13   USE in_out_manager  ! I/O manager
14   USE iom             ! for iom_put
15   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
16   USE timing          ! Timing
17   USE wrk_nemo        ! Wrk arrays
18   USE zdfmxl, ONLY : nmln
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC   zdf_convadj         ! called by step.F90
24   PUBLIC   zdf_convadj_init    ! called by nemogcm.f90
25
26   !
27   ! Namelist zdf_conadj
28   INTEGER nn_con                      ! Number of times to convect
29   CHARACTER(300) cn_convadj_mask      ! NetCDF mask file indicating points
30   LOGICAL , PUBLIC ::   ln_zdfconvadj ! convadj flag - used in step and tra_ldfiso
31   
32   REAL(wp) , PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: convadj_mask !: mask of points where adjustment applied
33   
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42
43CONTAINS
44
45   SUBROUTINE zdf_convadj( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE zdf_convadj  ***
48      !!                 
49      !! ** Purpose: Perform a convective adjustment of an unstable water column
50      !!               
51      !! ** Method: Given an unstable profile, the dense water will be moved
52      !!                down the water column to a depth at which it is stable, and the
53      !!                intervening boxes shuffled upwards, rather than using convection
54      !!                to stabilise the water column.
55      !!                Note that this requires a "cliff" to be dug out in bathymetry at
56      !!                overflow gridpoints.
57      !!                NB: This scheme modifies tsn in place to remove the
58      !!                    instabilities in these regions before other physics is
59      !!                    called.
60      !!
61      !! ** Action  :   T and S updated in static instability cases
62      !!----------------------------------------------------------------------
63      USE par_oce, ONLY: jp_tem, jp_sal, jpts
64
65      INTEGER, INTENT( in ) ::   kt   ! ocean time-step indexocean time step
66      INTEGER  ::   ji, jj, jk, jn, icon           ! dummy loop indices 
67      INTEGER  ::   itop, icona, iconb, ilevel
68      REAL(wp) , POINTER, DIMENSION(:) :: rhoup, rholo, trasum, tra
69      REAL(wp) ::   dztdif, tramix, rho, rho_insitu
70
71      IF( nn_timing == 1 )  CALL timing_start('zdf_convadj')
72      !
73      IF( kt == nit000 ) THEN
74         IF(lwp) WRITE(numout,*)
75         IF(lwp) WRITE(numout,*) 'zdf_convadj : Convective adjustment at overflow points'
76         IF(lwp) WRITE(numout,*) '~~~~~~~ '
77         IF(lwp) WRITE(numout,*)
78      ENDIF
79
80      !Allocation of arrays
81      CALL wrk_alloc( jpkm1, rhoup )
82      CALL wrk_alloc( jpkm1, rholo )
83      CALL wrk_alloc( jpts, tra )
84      CALL wrk_alloc( jpts, trasum )
85     
86      DO ji = 1,jpi !Do we need to include halos here?
87        DO jj = 1,jpj
88          IF (convadj_mask(ji,jj) .EQ. 1) THEN
89            DO icon = 1,nn_con ! Loop to repeat action nn_con times
90              rhoup(:)=0 
91              rholo(:)=0
92              itop=nmln(ji,jj)+1 ! Only apply below mixed layer
93              DO WHILE (itop .lt. mbkt(ji,jj))
94                ! Loop over depth and calculate potential density at W points (not surface)
95                ! *Do we need to worry about the before density or just the now?
96                DO jk = nmln(ji,jj),mbkt(ji,jj)-1
97                  rhoup(jk) = sigmai( tsn(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_sal), fsdepw(ji,jj,jk) )
98                  rholo(jk) = sigmai( tsn(ji,jj,jk+1,jp_tem), tsn(ji,jj,jk+1,jp_sal), fsdepw(ji,jj,jk) )
99                END DO
100 
101                ! Loop over depth from bottom to find shallowest
102                ! depth where column is unstable
103                icona = 0
104                DO jk = mbkt(ji,jj),itop,-1
105                  IF ( rhoup(jk) > rholo(jk) ) THEN
106                    icona = jk
107                    iconb = jk+1
108                  ENDIF
109                END DO
110
111                ! Loop over depth to find depth where this water is stable.
112                ! *Need to deal with key_vvl here.
113                IF (icona .NE. 0) THEN
114                  jk=iconb+1
115                  rho=1
116                  rho_insitu=0
117                  ilevel=iconb
118                  DO WHILE ((jk .LE. mbkt(ji,jj)) .AND. (rho .GT. rho_insitu))
119                    ilevel=jk-1
120                    rho = sigmai( tsn(ji,jj,icona,jp_tem), tsn(ji,jj,icona,jp_sal), fsdept(ji,jj,jk) )
121                    rho_insitu = sigmai( tsn(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_sal), fsdept(ji,jj,jk) )
122                    jk = jk+1
123                  END DO
124
125                  !Found the right level so replace water in lower box
126                  dztdif = fse3t(ji,jj,ilevel)-fse3t(ji,jj,icona)
127                  DO jn = 1, jpts ! Loop over tracers
128                    tra(jn) = tsn(ji,jj,ilevel,jn)
129                    trasum(jn) = tsn(ji,jj,icona,jn)*fse3t(ji,jj,icona) + tsn(ji,jj,ilevel,jn)*dztdif
130                    tramix = trasum(jn) / fse3t(ji,jj,ilevel)
131                    tsn(ji,jj,ilevel,jn) = tramix
132                  END DO
133
134
135                  !Shuffle up the intervening boxes
136                  DO jn = 1,jpts
137                    DO jk = ilevel-1,icona,-1
138                      dztdif = fse3t(ji,jj,jk)-fse3t(ji,jj,icona)
139                      trasum(jn) = tra(jn) * fse3t(ji,jj,icona) + tsn(ji,jj,jk,jn)*dztdif
140                      tramix = trasum(jn)/fse3t(ji,jj,jk)
141                      tra(jn) = tsn(ji,jj,jk,jn)
142                      tsn(ji,jj,jk,jn) = tramix
143                    END DO
144                  END DO
145
146                  itop = icona+1
147                ELSE
148                  itop=jpk
149                ENDIF
150              END DO ! End of WHILE loop
151            END DO !End loop over nn_con
152          END IF
153        END DO
154      END DO
155
156
157      CALL wrk_dealloc(jpkm1, rhoup, rholo)
158      CALL wrk_dealloc(jpts, tra, trasum)
159   END SUBROUTINE zdf_convadj
160
161
162   SUBROUTINE zdf_convadj_init
163      !!----------------------------------------------------------------------
164      !!                  ***  ROUTINE zdf_convadj  ***
165      !! ** Purpose: Initialise zdf_convadj module
166      !!
167      !! ** Method: Read namelist options
168      !!            Read mask from file
169      !!----------------------------------------------------------------------
170      USE iom
171
172      INTEGER               ::   imask
173
174      NAMELIST/namzdf_convadj/ nn_con, cn_convadj_mask, ln_convadj
175
176      REWIND ( numnam_ref )               !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
177      READ   ( numnam_ref, namzdf_convadj )
178
179
180      REWIND ( numnam_cfg )               !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
181      READ   ( numnam_cfg, namzdf_convadj )
182
183      IF(lwp) THEN                    !* Control print
184         WRITE(numout,*)
185         WRITE(numout,*) 'zdf_convadj_init : Convective adjustment scheme - initialisation'
186         WRITE(numout,*) '~~~~~~~~~~~~'
187         WRITE(numout,*) '   Namelist namzdf_convadj : set convadj parameters'
188         WRITE(numout,*) '      Use convective adjustment (cliff) scheme    ln_convadj= ', ln_convadj
189         WRITE(numout,*) '      Number of times to convect                  nn_con    = ', nn_con
190         WRITE(numout,*) '      NetCDF mask file of points to convect       cn_convadj_mask  = ',cn_convadj_mask
191      ENDIF
192
193      !Read in mask
194      IF( ln_convadj) THEN
195         ALLOCATE( convadj_mask(jpi,jpj) )
196         CALL iom_open ( cn_convadj_mask, imask )
197         CALL iom_get  ( imask, jpdom_autoglo, 'convadj_mask', convadj_mask)
198         CALL iom_close( imask )
199      ENDIF
200
201   END SUBROUTINE zdf_convadj_init
202
203   FUNCTION sigmai ( ptem, psal, pref )
204       !! --------------------------------------------------------------------
205       !! ** Purpose :   Compute the  density referenced to pref (ratio rho/rau0)
206       !!       from potential temperature and
207       !!      salinity fields using an equation of state defined through the
208       !!     namelist parameter neos.
209       !!
210       !! ** Method  :
211       !!       Jackett and McDougall (1994) equation of state.
212       !!         the in situ density is computed directly as a function of
213       !!         potential temperature relative to the surface (the opa t
214       !!         variable), salt and pressure (assuming no pressure variation
215       !!         along geopotential surfaces, i.e. the pressure p in decibars
216       !!         is approximated by the depth in meters.
217       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
218       !!              rhop(t,s)  = rho(t,s,0)
219       !!         with pressure                      p        decibars
220       !!              potential temperature         t        deg celsius
221       !!              salinity                      s        psu
222       !!              reference volumic mass        rau0     kg/m**3
223       !!              in situ volumic mass          rho      kg/m**3
224       !!              in situ density anomalie      prd      no units
225       !! --------------------------------------------------------------------
226 
227       !! * Arguments
228       REAL(wp), INTENT(in) :: ptem, psal
229       REAL(wp), INTENT(in) :: pref      !: reference pressure (meters or db)
230       REAL(wp)             :: sigmai
231 
232       REAL(kind=8),PARAMETER :: dpr4=4.8314d-4,dpd=-2.042967d-2 , dprau0 = 1000;
233
234       !! * local variables
235       INTEGER :: ji,jj 
236       REAL(wp) ::  dlrs
237       REAL(wp) ::  dlt, dls     
238       REAL(wp) :: dla,dla1,dlaw,dlb,dlb1,dlbw,dlc,dle,dlk0,dlkw 
239       REAL(wp) :: dlrhop, dlr1,dlr2,dlr3, dlref
240     
241       dlref = pref
242       sigmai = 0.d0
243       dlrs = SQRT( ABS( psal ) )
244
245! Convert T and S to double precision.
246        dlt=DBLE(ptem)
247        dls=DBLE(psal)
248
249
250! Compute the volumic mass of pure water at atmospheric pressure.
251        dlr1=((((6.536332d-9*dlt-1.120083d-6)&
252                *dlt+1.001685d-4)&
253               *dlt-9.095290d-3)&
254              *dlt+6.793952d-2)&
255             *dlt+999.842594d0
256
257! Compute the seawater volumic mass at atmospheric pressure.
258        dlr2=(((5.3875d-9*dlt-8.2467d-7)&
259               *dlt+7.6438d-5)&
260              *dlt-4.0899d-3)&
261             *dlt+0.824493d0
262
263        dlr3=(-1.6546d-6*dlt+1.0227d-4)&
264             *dlt-5.72466d-3
265
266! Compute the potential volumic mass (referenced to the surface).
267        dlrhop=(dpr4*dls+dlr3*dlrs+dlr2)*dls+dlr1
268
269! Compute the compression terms.
270        dle=(-3.508914d-8*dlt-1.248266d-8)&
271            *dlt-2.595994d-6
272
273        dlbw=(1.296821d-6*dlt-5.782165d-9)&
274             *dlt+1.045941d-4
275
276        dlb=dlbw+dle*dls
277
278        dlc=(-7.267926d-5*dlt+2.598241d-3 )&
279            *dlt+0.1571896d0
280
281        dlaw=((5.939910d-6*dlt+2.512549d-3)&
282              *dlt-0.1028859d0)&
283             *dlt-4.721788d0
284
285        dla=(dpd*dlrs+dlc)*dls+dlaw
286
287        dlb1=(-0.1909078d0*dlt+7.390729d0)&
288             *dlt-55.87545d0
289
290        dla1=((2.326469d-3*dlt+1.553190d0)&
291              *dlt-65.00517d0)&
292             *dlt+1044.077d0
293
294        dlkw=(((-1.361629d-4*dlt-1.852732d-2)&
295               *dlt-30.41638d0)&
296              *dlt+2098.925d0)&
297             *dlt+190925.6d0
298
299        dlk0=(dlb1*dlrs+dla1)*dls+dlkw
300
301! Compute the potential density anomaly.
302        sigmai=dlrhop/(1.0d0-dlref/(dlk0-dlref*(dla-dlref*dlb)))&
303                       -dprau0
304 
305    END FUNCTION sigmai
306
307END MODULE zdfconvadj
308
Note: See TracBrowser for help on using the repository browser.