source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90 @ 3421

Last change on this file since 3421 was 3421, checked in by charris, 9 years ago

#936 Changes as suggested by Gurvan et al (testing details in the ticket). Note that for ORCA1 results will change with either nn_closea=0 or 1 because the Caspian Sea is now coded as a closed sea.

  • Property svn:keywords set to Id
File size: 20.2 KB
Line 
1MODULE closea
2   !!======================================================================
3   !!                       ***  MODULE  closea  ***
4   !! Closed Seas  : specific treatments associated with closed seas
5   !!======================================================================
6   !! History :   8.2  !  00-05  (O. Marti)  Original code
7   !!             8.5  !  02-06  (E. Durand, G. Madec)  F90
8   !!             9.0  !  06-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
9   !!        NEMO 3.4  !  03-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   dom_clo    : modification of the ocean domain for closed seas cases
14   !!   sbc_clo    : Special handling of closed seas
15   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf)
16   !!   clo_ups    : set mixed centered/upstream scheme in closed sea (see traadv_cen2)
17   !!   clo_bat    : set to zero a field over closed sea (see domzrg)
18   !!----------------------------------------------------------------------
19   USE oce             ! dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
21   USE in_out_manager  ! I/O manager
22   USE sbc_oce         ! ocean surface boundary conditions
23   USE lib_fortran,    ONLY: glob_sum, DDPDD
24   USE lbclnk          ! lateral boundary condition - MPP exchanges
25   USE lib_mpp         ! MPP library
26   USE timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC dom_clo      ! routine called by domain module
32   PUBLIC sbc_clo      ! routine called by step module
33   PUBLIC clo_rnf      ! routine called by sbcrnf module
34   PUBLIC clo_ups      ! routine called in traadv_cen2(_jki) module
35   PUBLIC clo_bat      ! routine called in domzgr module
36
37   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea
38   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
39   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
40   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j)
41   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours
42   INTEGER, PUBLIC, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff
43
44   REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE dom_clo
56      !!---------------------------------------------------------------------
57      !!                  ***  ROUTINE dom_clo  ***
58      !!       
59      !! ** Purpose :   Closed sea domain initialization
60      !!
61      !! ** Method  :   if a closed sea is located only in a model grid point
62      !!                just the thermodynamic processes are applied.
63      !!
64      !! ** Action  :   ncsi1(), ncsj1() : south-west closed sea limits (i,j)
65      !!                ncsi2(), ncsj2() : north-east Closed sea limits (i,j)
66      !!                ncsir(), ncsjr() : Location of runoff
67      !!                ncsnr            : number of point where run-off pours
68      !!                ncstt            : Type of closed sea
69      !!                                   =0 spread over the world ocean
70      !!                                   =2 put at location runoff
71      !!----------------------------------------------------------------------
72      INTEGER ::   jc            ! dummy loop indices
73      !!----------------------------------------------------------------------
74     
75      IF(lwp) WRITE(numout,*)
76      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
77      IF(lwp) WRITE(numout,*)'~~~~~~~'
78
79      ! initial values
80      ncsnr(:) = 1  ;  ncsi1(:) = 1  ;  ncsi2(:) = 1  ;  ncsir(:,:) = 1
81      ncstt(:) = 0  ;  ncsj1(:) = 1  ;  ncsj2(:) = 1  ;  ncsjr(:,:) = 1
82
83      ! set the closed seas (in data domain indices)
84      ! -------------------
85
86      IF( cp_cfg == "orca" ) THEN
87         !
88         SELECT CASE ( jp_cfg )
89         !                                           ! =======================
90         CASE ( 1 )                                  ! ORCA_R1 configuration
91            !                                        ! =======================
92            ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea
93            ncsi1(1)   = 332  ; ncsj1(1)   = 203
94            ncsi2(1)   = 344  ; ncsj2(1)   = 235
95            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
96            !                                       
97            !                                        ! =======================
98         CASE ( 2 )                                  !  ORCA_R2 configuration
99            !                                        ! =======================
100            !                                            ! Caspian Sea
101            ncsnr(1)   =   1  ;  ncstt(1)   =   0           ! spread over the globe
102            ncsi1(1)   =  11  ;  ncsj1(1)   = 103
103            ncsi2(1)   =  17  ;  ncsj2(1)   = 112
104            ncsir(1,1) =   1  ;  ncsjr(1,1) =   1 
105            !                                            ! Great North American Lakes
106            ncsnr(2)   =   1  ;  ncstt(2)   =   2           ! put at St Laurent mouth
107            ncsi1(2)   =  97  ;  ncsj1(2)   = 107
108            ncsi2(2)   = 103  ;  ncsj2(2)   = 111
109            ncsir(2,1) = 110  ;  ncsjr(2,1) = 111
110            !                                            ! Black Sea 1 : west part of the Black Sea
111            ncsnr(3)   = 1    ; ncstt(3)   =   2            !            (ie west of the cyclic b.c.)
112            ncsi1(3)   = 174  ; ncsj1(3)   = 107            ! put in Med Sea
113            ncsi2(3)   = 181  ; ncsj2(3)   = 112
114            ncsir(3,1) = 171  ; ncsjr(3,1) = 106 
115            !                                            ! Black Sea 2 : est part of the Black Sea
116            ncsnr(4)   =   1  ;  ncstt(4)   =   2           !               (ie est of the cyclic b.c.)
117            ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! put in Med Sea
118            ncsi2(4)   =   6  ;  ncsj2(4)   = 112
119            ncsir(4,1) = 171  ;  ncsjr(4,1) = 106 
120            !                                        ! =======================
121         CASE ( 4 )                                  !  ORCA_R4 configuration
122            !                                        ! =======================
123            !                                            ! Caspian Sea
124            ncsnr(1)   =  1  ;  ncstt(1)   =  0 
125            ncsi1(1)   =  4  ;  ncsj1(1)   = 53 
126            ncsi2(1)   =  4  ;  ncsj2(1)   = 56
127            ncsir(1,1) =  1  ;  ncsjr(1,1) =  1
128            !                                            ! Great North American Lakes
129            ncsnr(2)   =  1  ;  ncstt(2)   =  2 
130            ncsi1(2)   = 49  ;  ncsj1(2)   = 55
131            ncsi2(2)   = 51  ;  ncsj2(2)   = 56
132            ncsir(2,1) = 57  ;  ncsjr(2,1) = 55
133            !                                            ! Black Sea
134            ncsnr(3)   =  4  ;  ncstt(3)   =  2 
135            ncsi1(3)   = 88  ;  ncsj1(3)   = 55 
136            ncsi2(3)   = 91  ;  ncsj2(3)   = 56
137            ncsir(3,1) = 86  ;  ncsjr(3,1) = 53
138            ncsir(3,2) = 87  ;  ncsjr(3,2) = 53 
139            ncsir(3,3) = 86  ;  ncsjr(3,3) = 52 
140            ncsir(3,4) = 87  ;  ncsjr(3,4) = 52
141            !                                            ! Baltic Sea
142            ncsnr(4)   =  1  ;  ncstt(4)   =  2
143            ncsi1(4)   = 75  ;  ncsj1(4)   = 59
144            ncsi2(4)   = 76  ;  ncsj2(4)   = 61
145            ncsir(4,1) = 84  ;  ncsjr(4,1) = 59 
146            !                                        ! =======================
147         CASE ( 025 )                                ! ORCA_R025 configuration
148            !                                        ! =======================
149            ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian + Aral sea
150            ncsi1(1)   = 1330 ; ncsj1(1)   = 645
151            ncsi2(1)   = 1400 ; ncsj2(1)   = 795
152            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
153            !                                       
154            ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Azov Sea
155            ncsi1(2)   = 1284 ; ncsj1(2)   = 722
156            ncsi2(2)   = 1304 ; ncsj2(2)   = 747
157            ncsir(2,1) = 1    ; ncsjr(2,1) = 1
158            !
159         END SELECT
160         !
161      ENDIF
162
163      ! convert the position in local domain indices
164      ! --------------------------------------------
165      DO jc = 1, jpncs
166         ncsi1(jc)   = mi0( ncsi1(jc) )
167         ncsj1(jc)   = mj0( ncsj1(jc) )
168
169         ncsi2(jc)   = mi1( ncsi2(jc) )   
170         ncsj2(jc)   = mj1( ncsj2(jc) ) 
171      END DO
172      !
173   END SUBROUTINE dom_clo
174
175
176   SUBROUTINE sbc_clo( kt )
177      !!---------------------------------------------------------------------
178      !!                  ***  ROUTINE sbc_clo  ***
179      !!                   
180      !! ** Purpose :   Special handling of closed seas
181      !!
182      !! ** Method  :   Water flux is forced to zero over closed sea
183      !!      Excess is shared between remaining ocean, or
184      !!      put as run-off in open ocean.
185      !!
186      !! ** Action  :   emp, emps   updated surface freshwater fluxes at kt
187      !!----------------------------------------------------------------------
188      INTEGER, INTENT(in) ::   kt   ! ocean model time step
189      !
190      INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices
191      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon
192      REAL(wp)            ::   zze2, ztmp, zcorr     !
193      COMPLEX(wp)         ::   ctmp 
194      REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace
195      !!----------------------------------------------------------------------
196      !
197      IF( nn_timing == 1 )  CALL timing_start('sbc_clo')
198      !                                                   !------------------!
199      IF( kt == nit000 ) THEN                             !  Initialisation  !
200         !                                                !------------------!
201         IF(lwp) WRITE(numout,*)
202         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
203         IF(lwp) WRITE(numout,*)'~~~~~~~'
204
205         surf(:) = 0.e0_wp
206         !
207         surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean
208         !
209         !                                        ! surface of closed seas
210         IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
211            DO jc = 1, jpncs
212               ctmp = CMPLX( 0.e0, 0.e0, wp )
213               DO jj = ncsj1(jc), ncsj2(jc)
214                  DO ji = ncsi1(jc), ncsi2(jc)
215                     ztmp = e1e2t(ji,jj) * tmask_i(ji,jj)
216                     CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
217                  END DO
218               END DO
219               IF( lk_mpp )   CALL mpp_sum( ctmp )
220               surf(jc) = REAL(ctmp,wp)
221            END DO
222         ELSE                                          ! Standard calculation           
223            DO jc = 1, jpncs
224               DO jj = ncsj1(jc), ncsj2(jc)
225                  DO ji = ncsi1(jc), ncsi2(jc)
226                     surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas
227                  END DO
228               END DO
229            END DO
230            IF( lk_mpp )   CALL mpp_sum ( surf, jpncs+1 )       ! mpp: sum over all the global domain
231         ENDIF
232
233         IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
234         DO jc = 1, jpncs
235            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
236         END DO
237
238         ! jpncs+1 : surface of sea, closed seas excluded
239         DO jc = 1, jpncs
240            surf(jpncs+1) = surf(jpncs+1) - surf(jc)
241         END DO           
242         !
243      ENDIF
244      !                                                   !--------------------!
245      !                                                   !  update emp, emps  !
246      zfwf = 0.e0_wp                                      !--------------------!
247      IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
248         DO jc = 1, jpncs
249            ctmp = CMPLX( 0.e0, 0.e0, wp )
250            DO jj = ncsj1(jc), ncsj2(jc)
251               DO ji = ncsi1(jc), ncsi2(jc)
252                  ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)
253                  CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
254               END DO 
255            END DO
256            IF( lk_mpp )   CALL mpp_sum( ctmp )
257            zfwf(jc) = REAL(ctmp,wp)
258         END DO
259      ELSE                                          ! Standard calculation           
260         DO jc = 1, jpncs
261            DO jj = ncsj1(jc), ncsj2(jc)
262               DO ji = ncsi1(jc), ncsi2(jc)
263                  zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 
264               END DO 
265            END DO
266         END DO
267         IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain
268      ENDIF
269
270      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration
271         zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp
272         zfwf(3) = zze2
273         zfwf(4) = zze2
274      ENDIF
275
276      zcorr = 0._wp
277
278      DO jc = 1, jpncs
279         !
280         ! The following if avoids the redistribution of the round off
281         IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN
282            !
283            IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean
284               emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1)
285               emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1)
286               ! accumulate closed seas correction
287               zcorr     = zcorr     + zfwf(jc) / surf(jpncs+1)
288               !
289            ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared
290               IF ( zfwf(jc) <= 0.e0_wp ) THEN
291                   DO jn = 1, ncsnr(jc)
292                     ji = mi0(ncsir(jc,jn))
293                     jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
294                     IF (      ji > 1 .AND. ji < jpi   &
295                         .AND. jj > 1 .AND. jj < jpj ) THEN
296                         emp (ji,jj) = emp (ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )
297                         emps(ji,jj) = emps(ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )
298                     ENDIF
299                   END DO
300               ELSE
301                   emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1)
302                   emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1)
303                   ! accumulate closed seas correction
304                   zcorr     = zcorr     + zfwf(jc) / surf(jpncs+1)
305               ENDIF
306            ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location
307               DO jn = 1, ncsnr(jc)
308                  ji = mi0(ncsir(jc,jn))
309                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
310                  IF(      ji > 1 .AND. ji < jpi    &
311                     .AND. jj > 1 .AND. jj < jpj ) THEN
312                     emp (ji,jj) = emp (ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) )
313                     emps(ji,jj) = emps(ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) )
314                  ENDIF
315               END DO
316            ENDIF 
317            !
318            DO jj = ncsj1(jc), ncsj2(jc)
319               DO ji = ncsi1(jc), ncsi2(jc)
320                  emp (ji,jj) = emp (ji,jj) - zfwf(jc) / surf(jc)
321                  emps(ji,jj) = emps(ji,jj) - zfwf(jc) / surf(jc)
322               END DO 
323            END DO 
324            !
325         END IF
326      END DO
327
328      IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas
329         DO jc = 1, jpncs                  ! only if it is large enough
330            DO jj = ncsj1(jc), ncsj2(jc)
331               DO ji = ncsi1(jc), ncsi2(jc)
332                  emp (ji,jj) = emp (ji,jj) - zcorr
333                  emps(ji,jj) = emps(ji,jj) - zcorr
334               END DO 
335             END DO
336          END DO
337      ENDIF
338      !
339      emp (:,:) = emp (:,:) * tmask(:,:,1)
340      emps(:,:) = emps(:,:) * tmask(:,:,1)
341      !
342      CALL lbc_lnk( emp , 'T', 1._wp )
343      CALL lbc_lnk( emps, 'T', 1._wp )
344      !
345      IF( nn_timing == 1 )  CALL timing_stop('sbc_clo')
346      !
347   END SUBROUTINE sbc_clo
348
349
350   SUBROUTINE clo_rnf( p_rnfmsk )
351      !!---------------------------------------------------------------------
352      !!                  ***  ROUTINE sbc_rnf  ***
353      !!                   
354      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
355      !!                to be the same as river mouth grid-points
356      !!
357      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
358      !!                at the closed sea outflow grid-point.
359      !!
360      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
361      !!----------------------------------------------------------------------
362      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
363      !
364      INTEGER  ::   jc, jn      ! dummy loop indices
365      INTEGER  ::   ii, ij      ! temporary integer
366      !!----------------------------------------------------------------------
367      !
368      DO jc = 1, jpncs
369         IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows
370             DO jn = 1, 4
371               ii = mi0( ncsir(jc,jn) )
372               ij = mj0( ncsjr(jc,jn) )
373               p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0_wp )
374            END DO
375         ENDIF
376      END DO 
377      !
378   END SUBROUTINE clo_rnf
379
380   
381   SUBROUTINE clo_ups( p_upsmsk )
382      !!---------------------------------------------------------------------
383      !!                  ***  ROUTINE sbc_rnf  ***
384      !!                   
385      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
386      !!                to be the same as river mouth grid-points
387      !!
388      !! ** Method  :   set to 0.5 the upstream mask (upsmsk, see traadv_cen2
389      !!                module) over the closed seas.
390      !!
391      !! ** Action  :   update (p_)upsmsk (set 0.5 over closed seas)
392      !!----------------------------------------------------------------------
393      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_upsmsk   ! upstream mask (upsmsk array)
394      !
395      INTEGER  ::   jc, ji, jj      ! dummy loop indices
396      !!----------------------------------------------------------------------
397      !
398      DO jc = 1, jpncs
399         DO jj = ncsj1(jc), ncsj2(jc)
400            DO ji = ncsi1(jc), ncsi2(jc)
401               p_upsmsk(ji,jj) = 0.5_wp         ! mixed upstream/centered scheme over closed seas
402            END DO
403         END DO
404       END DO 
405       !
406   END SUBROUTINE clo_ups
407   
408     
409   SUBROUTINE clo_bat( pbat, kbat )
410      !!---------------------------------------------------------------------
411      !!                  ***  ROUTINE clo_bat  ***
412      !!                   
413      !! ** Purpose :   suppress closed sea from the domain
414      !!
415      !! ** Method  :   set to 0 the meter and level bathymetry (given in
416      !!                arguments) over the closed seas.
417      !!
418      !! ** Action  :   set pbat=0 and kbat=0 over closed seas
419      !!----------------------------------------------------------------------
420      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pbat   ! bathymetry in meters (bathy array)
421      INTEGER , DIMENSION(jpi,jpj), INTENT(inout) ::   kbat   ! bathymetry in levels (mbathy array)
422      !
423      INTEGER  ::   jc, ji, jj      ! dummy loop indices
424      !!----------------------------------------------------------------------
425      !
426      DO jc = 1, jpncs
427         DO jj = ncsj1(jc), ncsj2(jc)
428            DO ji = ncsi1(jc), ncsi2(jc)
429               pbat(ji,jj) = 0._wp   
430               kbat(ji,jj) = 0   
431            END DO
432         END DO
433       END DO 
434       !
435   END SUBROUTINE clo_bat
436
437   !!======================================================================
438END MODULE closea
439
Note: See TracBrowser for help on using the repository browser.