source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90 @ 3632

Last change on this file since 3632 was 3632, checked in by acc, 9 years ago

Branch dev_NOC_2012_r3555. #1006. Step 9: Merge in trunk changes between revision 3385 and 3452

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