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.
closea.F90 in branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_r8600_closea_rewrite/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90 @ 9017

Last change on this file since 9017 was 9017, checked in by davestorkey, 6 years ago

UKMO/dev_r8600_closea_rewrite branch : commit code

File size: 19.0 KB
Line 
1MODULE closea
2   !!======================================================================
3   !!                   ***  MODULE  usrdef_closea  ***
4   !!
5   !! User define : specific treatments associated with closed seas
6   !!======================================================================
7   !! History :   8.2  !  2000-05  (O. Marti)  Original code
8   !!   NEMO      1.0  !  2002-06  (E. Durand, G. Madec)  F90
9   !!             3.0  !  2006-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
10   !!             3.4  !  2014-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility
11   !!             4.0  !  2016-06  (G. Madec)  move to usrdef_closea, remove clo_ups
12   !!             4.0  !  2017-12  (D. Storkey) new formulation based on masks read from file
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   dom_clo    : modification of the ocean domain for closed seas cases
17   !!   sbc_clo    : Special handling of closed seas
18   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf)
19   !!   clo_bat    : set to zero a field over closed sea (see domzrg)
20   !!----------------------------------------------------------------------
21   USE oce             ! dynamics and tracers
22   USE dom_oce         ! ocean space and time domain
23   USE phycst          ! physical constants
24   USE sbc_oce         ! ocean surface boundary conditions
25   USE iom             ! I/O routines
26   !
27   USE in_out_manager  ! I/O manager
28   USE lib_fortran,    ONLY: glob_sum
29   USE lbclnk          ! lateral boundary condition - MPP exchanges
30   USE lib_mpp         ! MPP library
31   USE timing
32   USE wrk_nemo        ! Memory allocation
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC dom_clo      ! called by domain module
38   PUBLIC sbc_clo      ! called by sbcmod module
39   PUBLIC clo_rnf      ! called by sbcrnf module
40   PUBLIC clo_bat      ! called in domain module
41
42   INTEGER, PUBLIC :: jncs !: number of closed seas (inferred from closea_mask field)
43   INTEGER, PUBLIC :: jncsr !: number of closed seas rnf mappings (inferred from closea_mask_rnf field)
44   INTEGER, PUBLIC :: jncse !: number of closed seas empmr mappings (inferred from closea_mask_empmr field)
45   
46   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask       !: mask of integers defining closed seas
47   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask_rnf   !: mask of integers defining closed seas rnf mappings
48   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::  closea_mask_empmr !: mask of integers defining closed seas empmr mappings
49   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surf         !: closed sea surface areas
50                                                                  !: (and residual global surface area)
51   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surfr        !: closed sea target rnf surface areas
52   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:)  ::   surfe        !: closed sea target empmr surface areas
53
54   !! * Substitutions
55#  include "vectopt_loop_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE dom_clo( k_bot )
64      !!---------------------------------------------------------------------
65      !!                  ***  ROUTINE dom_clo  ***
66      !!       
67      !! ** Purpose :   Closed sea domain initialization
68      !!
69      !! ** Method  :   if a closed sea is located only in a model grid point
70      !!                just the thermodynamic processes are applied.
71      !!
72      !! ** Action  :   Read closea mask fields from closea_masks.nc and infer
73      !!                number of closed seas from closea_mask field.
74      !!                closea_mask       : integer values defining closed seas (or groups of closed seas)
75      !!                closea_mask_rnf   : integer values defining mappings from closed seas or groups of
76      !!                                    closed seas to a runoff area for downwards flux only.
77      !!                closea_mask_empmr : integer values defining mappings from closed seas or groups of
78      !!                                    closed seas to a runoff area for net fluxes.
79      !!
80      !!                Python code to generate the closea_masks.nc file from the old-style indices
81      !!                definitions is available at TOOLS/DOMAINcfg/make_closea_masks.py
82      !!----------------------------------------------------------------------
83      INTEGER, DIMENSION(:,:), INTENT(in)  ::   k_bot   ! ocean last level index (for masking input fields)
84      INTEGER ::   inum    ! input file identifier
85      INTEGER ::   ierr    ! error code
86      INTEGER ::   id      ! netcdf variable ID
87      REAL(wp), POINTER, DIMENSION(:,:) :: zdata_in ! temporary real array for input
88      !!----------------------------------------------------------------------
89      !
90      IF(lwp) WRITE(numout,*)
91      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
92      IF(lwp) WRITE(numout,*)'~~~~~~~'
93      !
94      CALL wrk_alloc( jpi,jpj, zdata_in)
95      !
96      ! read the closed seas masks
97      ! --------------------------
98      !
99      ALLOCATE( closea_mask(jpi,jpj) , STAT=ierr )
100      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask array')
101      ALLOCATE( closea_mask_rnf(jpi,jpj) , STAT=ierr )
102      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_rnf array')
103      ALLOCATE( closea_mask_empmr(jpi,jpj) , STAT=ierr )
104      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dom_clo: failed to allocate closea_mask_empmr array')
105
106      closea_mask(:,:) = 0
107      closea_mask_rnf(:,:) = 0
108      closea_mask_empmr(:,:) = 0
109
110      CALL iom_open( 'closea_masks.nc', inum )
111
112      zdata_in(:,:) = 0.0
113      CALL iom_get ( inum, jpdom_data, 'closea_mask', zdata_in )
114      closea_mask(:,:) = NINT(zdata_in(:,:))
115      ! necessary because fill values in input fields can confuse things
116      ! we can't multiply by tmask here because it isn't defined yet
117      WHERE( k_bot(:,:) == 0 ) closea_mask(:,:) = 0
118
119      id = iom_varid(inum, 'closea_mask_rnf', ldstop = .false.)
120      IF(lwp) WRITE(numout,*) 'closea_mask_rnf, id : ',id
121      IF( id > 0 ) THEN
122         CALL iom_get ( inum, jpdom_data, 'closea_mask_rnf', zdata_in )
123         closea_mask_rnf(:,:) = NINT(zdata_in(:,:))
124         WHERE( k_bot(:,:) == 0 ) closea_mask_rnf(:,:) = 0
125      ENDIF
126 
127      id = iom_varid(inum, 'closea_mask_empmr', ldstop = .false.)
128      IF( id > 0 ) THEN
129         CALL iom_get ( inum, jpdom_data, 'closea_mask_empmr', zdata_in )
130         closea_mask_empmr(:,:) = NINT(zdata_in(:,:))
131         WHERE( k_bot(:,:) == 0 ) closea_mask_empmr(:,:) = 0
132      ENDIF
133
134      CALL iom_close( inum )
135
136      ! number of closed seas = global maximum value in closea_mask field
137      jncs = maxval(closea_mask(:,:))
138      IF( lk_mpp ) CALL mpp_max(jncs)
139      IF( lwp ) WRITE(numout,*) 'Number of closed seas : ',jncs
140      !
141      ! number of closed seas rnf mappings = global maximum in closea_mask_rnf field
142      jncsr = maxval(closea_mask_rnf(:,:))
143      IF( lk_mpp ) CALL mpp_max(jncsr)
144      IF( lwp ) WRITE(numout,*) 'Number of closed seas rnf mappings : ',jncsr
145      !
146      ! number of closed seas empmr mappings = global maximum value in closea_mask_empmr field
147      jncse = maxval(closea_mask_empmr(:,:))
148      IF( lk_mpp ) CALL mpp_max(jncse)
149      IF( lwp ) WRITE(numout,*) 'Number of closed seas empmr mappings : ',jncse
150      !
151      CALL wrk_dealloc(jpi,jpj, zdata_in)
152      !
153   END SUBROUTINE dom_clo
154
155
156   SUBROUTINE sbc_clo( kt )
157      !!---------------------------------------------------------------------
158      !!                  ***  ROUTINE sbc_clo  ***
159      !!                   
160      !! ** Purpose :   Special handling of closed seas
161      !!
162      !! ** Method  :   Water flux is forced to zero over closed sea
163      !!      Excess is shared between remaining ocean, or
164      !!      put as run-off in open ocean.
165      !!
166      !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt
167      !!----------------------------------------------------------------------
168      INTEGER         , INTENT(in   ) ::   kt       ! ocean model time step
169      !
170      INTEGER             ::   ierr
171      INTEGER             ::   jc, jcr, jce   ! dummy loop indices
172      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon
173      REAL(wp)            ::   zfwf_total, zcoef, zcoef1         !
174      REAL(wp), POINTER, DIMENSION(:)   ::   zfwf, zfwfr, zfwfe     ! freshwater fluxes over closed seas
175      REAL(wp), POINTER, DIMENSION(:,:) ::   ztmp2d   ! 2D workspace
176      !!----------------------------------------------------------------------
177      !
178      IF( nn_timing == 1 )  CALL timing_start('sbc_clo')
179      !
180      CALL wrk_alloc( jncs, zfwf )
181      ! jncsr and jncse can be zero so add 1 to avoid allocating zero-length arrays
182      CALL wrk_alloc( jncsr+1, zfwfr )
183      CALL wrk_alloc( jncse+1, zfwfe )
184      CALL wrk_alloc( jpi,jpj, ztmp2d )
185      !                                                   !------------------!
186      IF( kt == nit000 ) THEN                             !  Initialisation  !
187         !                                                !------------------!
188         IF(lwp) WRITE(numout,*)
189         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
190         IF(lwp) WRITE(numout,*)'~~~~~~~'
191
192         ALLOCATE( surf(jncs+1) , STAT=ierr )
193         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surf array')
194         surf(:) = 0.e0_wp
195         !
196         ! jncsr can be zero so add 1 to avoid allocating zero-length array
197         ALLOCATE( surfr(jncsr+1) , STAT=ierr )
198         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfr array')
199         surfr(:) = 0.e0_wp
200         !
201         ! jncse can be zero so add 1 to avoid allocating zero-length array
202         ALLOCATE( surfe(jncse+1) , STAT=ierr )
203         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sbc_clo: failed to allocate surfe array')
204         surfe(:) = 0.e0_wp
205         !
206         surf(jncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean
207         !
208         !                                        ! surface areas of closed seas
209         DO jc = 1, jncs
210            ztmp2d(:,:) = 0.e0_wp
211            WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:)
212            surf(jc) = glob_sum( ztmp2d(:,:) )
213         END DO
214         !
215         ! jncs+1 : surface area of global ocean, closed seas excluded
216         surf(jncs+1) = surf(jncs+1) - SUM(surf(1:jncs))
217         !
218         !                                        ! surface areas of rnf target areas
219         DO jcr = 1, jncsr
220            ztmp2d(:,:) = 0.e0_wp
221            WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:)
222            surfr(jcr) = glob_sum( ztmp2d(:,:) )
223         END DO
224         !
225         !                                        ! surface areas of empmr target areas
226         DO jce = 1, jncse
227            ztmp2d(:,:) = 0.e0_wp
228            WHERE( closea_mask_empmr(:,:) == jcr .and. closea_mask(:,:) == 0 ) ztmp2d(:,:) = e1e2t(:,:) * tmask_i(:,:)
229            surfe(jcr) = glob_sum( ztmp2d(:,:) )
230         END DO
231         !
232         IF(lwp) WRITE(numout,*)'     Closed sea surface areas (km2)'
233         DO jc = 1, jncs
234            IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jc, surf(jc) * 1.0e-6
235         END DO
236         IF(lwp) WRITE(numout,FMT='(A,ES12.2)') 'Global surface area excluding closed seas (km2): ', surf(jncs+1) * 1.0e-6
237         !
238         IF(jncsr > 0) THEN
239            IF(lwp) WRITE(numout,*)'     Closed sea target rnf surface areas (km2)'
240            DO jcr = 1, jncsr
241               IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jcr, surfr(jcr) * 1.0e-6
242            END DO
243         ENDIF
244         !
245         IF(jncse > 0) THEN
246            IF(lwp) WRITE(numout,*)'     Closed sea target empmr surface areas (km2)'
247            DO jce = 1, jncse
248               IF(lwp) WRITE(numout,FMT='(1I3,5X,ES12.2)') jce, surfe(jce) * 1.0e-6
249            END DO
250         ENDIF
251      ENDIF
252      !
253      !                                                      !--------------------!
254      !                                                      !  update emp        !
255      !                                                      !--------------------!
256
257      zfwf_total = 0._wp
258
259      !
260      ! 1. Work out total freshwater fluxes over closed seas from EMP - RNF.
261      !
262      zfwf(:) = 0.e0_wp           
263      DO jc = 1, jncs
264         ztmp2d(:,:) = 0.e0_wp
265         WHERE( closea_mask(:,:) == jc ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:)
266         zfwf(jc) = glob_sum( ztmp2d(:,:) )
267      END DO
268      zfwf_total = SUM(zfwf)
269
270      !
271      ! 2. Work out total FW fluxes over rnf source areas and add to rnf target areas. If jncsr is zero does not loop.
272      !    Where zfwf is negative add flux at specified runoff points and subtract from fluxes for global redistribution.
273      !    Where positive leave in global redistribution total.
274      !
275      zfwfr(:) = 0.e0_wp           
276      DO jcr = 1, jncsr
277         !
278         ztmp2d(:,:) = 0.e0_wp
279         WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:)
280         zfwfr(jcr) = glob_sum( ztmp2d(:,:) )
281         IF(lwp) WRITE(numout,*) 'closea runoff output: jcr, zfwfr(jcr), ABS(zfwfr(jcr)/surf(jncs+1) : ',jcr,zfwfr(jcr),ABS(zfwfr(jcr)/surf(jncs+1))
282         !
283         ! The following if avoids the redistribution of the round off
284         IF ( ABS(zfwfr(jcr) / surf(jncs+1) ) > rsmall) THEN
285            !
286            ! Add residuals to target runoff points if negative and subtract from total to be added globally
287            IF(lwp) WRITE(numout,*) 'closea runoff output: passed roundoff test!'
288            IF( zfwfr(jcr) < 0.0 ) THEN
289               zfwf_total = zfwf_total - zfwfr(jcr)
290               zcoef    = zfwfr(jcr) / surfr(jcr)
291               zcoef1   = rcp * zcoef
292               IF(lwp) WRITE(numout,*) 'closea runoff output: jcr, zcoef: ',jcr, zcoef
293               WHERE( closea_mask_rnf(:,:) == jcr .and. closea_mask(:,:) == 0.0)
294                  emp(:,:) = emp(:,:) - zcoef
295                  qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:)
296               ENDWHERE
297            ENDIF
298            !
299         ENDIF
300      END DO
301   
302      !
303      ! 3. Work out total fluxes over empmr source areas and add to empmr target areas. If jncse is zero does not loop.
304      !
305      zfwfe(:) = 0.e0_wp           
306      DO jce = 1, jncse
307         !
308         ztmp2d(:,:) = 0.e0_wp
309         WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) > 0 ) ztmp2d(:,:) = e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:)
310         zfwfe(jce) = glob_sum( ztmp2d(:,:) )
311         !
312         ! The following if avoids the redistribution of the round off
313         IF ( ABS(zfwfe(jce) / surf(jncs+1) ) > rsmall) THEN
314            !
315            ! Add residuals to runoff points and subtract from total to be added globally
316            zfwf_total = zfwf_total - zfwfe(jce)
317            zcoef    = zfwfe(jce) / surfe(jce)
318            zcoef1   = rcp * zcoef
319            WHERE( closea_mask_empmr(:,:) == jce .and. closea_mask(:,:) == 0.0)
320               emp(:,:) = emp(:,:) - zcoef
321               qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:)
322            ENDWHERE
323            !
324         ENDIF
325      END DO
326   
327      !
328      ! 4. Spread residual flux over global ocean.
329      !
330      ! The following if avoids the redistribution of the round off
331      IF ( ABS(zfwf_total / surf(jncs+1) ) > rsmall) THEN
332         zcoef    = zfwf_total / surf(jncs+1)
333         zcoef1   = rcp * zcoef
334         IF(lwp) WRITE(numout,*) 'closea global addition: zfwf_total, zcoef: ', zfwf_total, zcoef
335         WHERE( closea_mask(:,:) == 0 )
336            emp(:,:) = emp(:,:) + zcoef
337            qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
338         ENDWHERE
339      ENDIF
340
341      !
342      ! 5. Subtract area means from emp (and qns) over closed seas to give zero mean FW flux over each sea.
343      !
344      DO jc = 1, jncs
345         ! The following if avoids the redistribution of the round off
346         IF ( ABS(zfwf(jc) / surf(jncs+1) ) > rsmall) THEN
347            !
348            ! Subtract residuals from fluxes over closed sea
349            zcoef    = zfwf(jc) / surf(jc)
350            zcoef1   = rcp * zcoef
351            IF(lwp) WRITE(numout,*) 'closea subtract mean: jc, zfwf(jc), zcoef: ',jc, zfwf(jc), zcoef
352            WHERE( closea_mask(:,:) == jc )
353               emp(:,:) = emp(:,:) - zcoef
354               qns(:,:) = qns(:,:) + zcoef1 * sst_m(:,:)
355            ENDWHERE
356            !
357         ENDIF
358      END DO
359      !
360      emp (:,:) = emp (:,:) * tmask(:,:,1)
361      !
362      CALL lbc_lnk( emp , 'T', 1._wp )
363      !
364      CALL wrk_dealloc( jncs   , zfwf )
365      CALL wrk_dealloc( jncsr+1  , zfwfr )
366      CALL wrk_dealloc( jncse+1  , zfwfe )
367      CALL wrk_dealloc( jpi,jpj, ztmp2d )
368      !
369      IF( nn_timing == 1 )  CALL timing_stop('sbc_clo')
370      !
371   END SUBROUTINE sbc_clo
372
373   SUBROUTINE clo_rnf( p_rnfmsk )
374      !!---------------------------------------------------------------------
375      !!                  ***  ROUTINE sbc_rnf  ***
376      !!                   
377      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
378      !!                to be the same as river mouth grid-points
379      !!
380      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
381      !!                at the closed sea outflow grid-point.
382      !!
383      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
384      !!----------------------------------------------------------------------
385      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
386      !!----------------------------------------------------------------------
387      !
388      WHERE( ( closea_mask_rnf(:,:) > 0 .or. closea_mask_empmr(:,:) > 0 ) .and. closea_mask(:,:) == 0 )
389         p_rnfmsk(:,:) = MAX( p_rnfmsk(:,:), 1.0_wp )
390      ENDWHERE
391      !
392   END SUBROUTINE clo_rnf
393   
394     
395   SUBROUTINE clo_bat( k_top, k_bot )
396      !!---------------------------------------------------------------------
397      !!                  ***  ROUTINE clo_bat  ***
398      !!                   
399      !! ** Purpose :   suppress closed sea from the domain
400      !!
401      !! ** Method  :   set first and last ocean level to 0 over the closed seas.
402      !!
403      !! ** Action  :   set pbat=0 and kbat=0 over closed seas
404      !!----------------------------------------------------------------------
405      INTEGER, DIMENSION(:,:), INTENT(inout) ::   k_top, k_bot   ! ocean first and last level indices
406      !!----------------------------------------------------------------------
407      !
408      WHERE( closea_mask(:,:) > 0 )
409         k_top(:,:) = 0   
410         k_bot(:,:) = 0   
411      ENDWHERE
412      !
413   END SUBROUTINE clo_bat
414
415   !!======================================================================
416END MODULE closea
417
Note: See TracBrowser for help on using the repository browser.