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.
domclo.F90 in NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src – NEMO

source: NEMO/branches/2019/ENHANCE-02_ISF_domcfg/src/domclo.F90 @ 11568

Last change on this file since 11568 was 11201, checked in by mathiot, 5 years ago

ENHANCE-03_domcfg : add management of closed seas in domain cfg by flood filling and lat/lon seed instead of i/j box definition (ticket #2143)
File size: 14.2 KB
Line 
1MODULE domclo
2   !!======================================================================
3   !!                       ***  MODULE  domclo  ***
4   !! domclo : definition of closed sea mask needed by NEMO
5   !!======================================================================
6   !! History : NEMO 4.0 ! 07-2019 (P. Mathiot) original version
7   !!----------------------------------------------------------------------
8   !!----------------------------------------------------------------------
9   !!   dom_clo    : definition of closed sea mask needed by NEMO
10   !!----------------------------------------------------------------------
11   USE dom_oce         ! ocean space and time domain
12   USE domngb          ! closest point algorithm
13   USE phycst          ! rpi, rad, ra
14   USE domutil         ! flood filling algorithm (fill_pool)
15   USE in_out_manager  ! I/O manager
16   USE lbclnk          ! lateral boundary condition - MPP exchanges
17   USE lib_mpp
18   USE lib_fortran
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC dom_clo      ! routine called by domain module
24
25   !! * Substitutions
26#  include "vectopt_loop_substitute.h90"
27   !!----------------------------------------------------------------------
28   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
29   !! $Id: closea.F90 7213 2016-11-09 09:07:42Z timgraham $
30   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
31   !!----------------------------------------------------------------------
32CONTAINS
33
34   SUBROUTINE dom_clo
35      !!---------------------------------------------------------------------
36      !!                  ***  ROUTINE dom_clo  ***
37      !!       
38      !! ** Purpose :   Closed sea mask definition
39      !!
40      !! ** Method  :   A flood filling algorithm is used to detect lake and open ocean
41      !!                Namelist provide seed for each lake and condition on how to
42      !!                manage net evap/runoff.
43      !!
44      !! ** Action  :  - compute open ocean mask
45      !!               - loop over all the lake in namelist
46      !!               - define lake msk (msk_csglo/emp/rnf)
47      !!               - define river mouth mask for each lake or group of lake with the same river mouth
48      !!                 (msk_csgrpglo/emp/rnf)
49      !!               
50      !!----------------------------------------------------------------------
51
52      TYPE closea 
53         CHARACTER(256) :: cname                       ! name
54         REAL(wp)       :: rlonsrc                     ! seed src location (lon)
55         REAL(wp)       :: rlatsrc                     ! seed src location (lat)
56         REAL(wp)       :: rlontrg                     ! seed trg location (lon)
57         REAL(wp)       :: rlattrg                     ! seed trg location (lat)
58         CHARACTER(256) :: cloctrg                     ! where water is spread
59         CHARACTER(256) :: cschtrg                     ! how   water is spread
60         REAL(wp)       :: radtrg                      ! radius of closed sea river mouth
61         INTEGER        :: idtrg                       ! target id in case multiple lakes for the same river mouth
62      END TYPE
63         
64      INTEGER, PARAMETER :: jp_closea = 98 ! number maximal of closea
65      INTEGER :: ji, jj, jcs, jsch  ! loop indexes
66      INTEGER :: jglo, jemp, jrnf   ! closed sea indexes for global, emp, rnf spreading
67      INTEGER :: jiseed, jjseed     ! seed indexes
68      INTEGER :: jglo, jloc, jcoast ! counter
69      INTEGER :: ios
70      INTEGER :: nn_closea          ! number of closed seas
71
72      REAL(wp) :: zdistseed         ! distance to seed
73      REAL(wp) :: zarea             ! river mouth area
74      REAL(wp) :: rn_lon_opnsea, rn_lat_opnsea ! open sea seed
75      REAL(wp), DIMENSION(1)       :: zchk, zradtrg
76      REAL(wp), DIMENSION(jpi,jpj) :: zmsksrc, zmsktrg, zmsk_coastline, zdist
77
78      CHARACTER(256) :: csch, cloc        ! scheme name for water spreading (glo, rnf, emp)
79      TYPE(closea)  , DIMENSION(jp_closea)   :: sn_lake   ! lake properties
80
81      LOGICAL :: lskip     ! flag in case lake seed on land or already filled (...)
82
83      !!----------------------------------------------------------------------
84      !! 0 : Read namelist for closed sea definition
85      !!----------------------------------------------------------------------
86      IF(lwp) WRITE(numout,*)
87      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
88      IF(lwp) WRITE(numout,*)'~~~~~~~'
89
90      NAMELIST/namclo/ rn_lon_opnsea, rn_lat_opnsea, nn_closea, sn_lake
91      !!---------------------------------------------------------------------
92      PRINT *, rn_lon_opnsea, rn_lat_opnsea, nn_closea
93      !!
94      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
95      READ  ( numnam_ref, namclo, IOSTAT = ios, ERR = 901 )
96901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namclo in reference namelist', lwp )
97      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
98      READ  ( numnam_cfg, namclo, IOSTAT = ios, ERR = 902 )
99902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namclo in configuration namelist', lwp )
100      IF(lwm) WRITE ( numond, namclo )
101
102      !!----------------------------------------------------------------------
103      !! 1 : Define open ocean and closed sea masks
104      !!----------------------------------------------------------------------
105
106      !! 1.1 get closest point to the lat/lon seed
107      msk_opnsea(:,:) = ssmask(:,:)
108      CALL dom_ngb(rn_lon_opnsea, rn_lat_opnsea, jiseed, jjseed, zdistseed, 'T')
109
110      !! 1.2 fill select cell to -1
111      CALL fill_pool( jiseed, jjseed, msk_opnsea, -1._wp )
112
113      !! sanity check on the msk value (closea mask should be 99 otherwise, lake already processed)
114      !! check if some lake are not connected
115      zchk = 0._wp
116      IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) zchk = ssmask(mi0(jiseed),mj0(jjseed))
117      CALL mpp_max('domclo',zchk)
118      IF (zchk(1) == 0._wp) CALL ctl_stop( 'STOP', 'open sea seed is on land, please update namelist (rn_lon_opnsea,rn_lat_opnsea)' ) 
119
120      !! 1.3 set to 0 everything >0 and revert mask
121      IF (lwp) THEN
122         WRITE(numout,*)
123         WRITE(numout,*)'open ocean: '
124         WRITE(numout,'(a,2f7.2)')'    lon/lat seed to detect main ocean is: ', rn_lon_opnsea, rn_lat_opnsea
125         WRITE(numout,'(a,2i7)'  )'    i/j     seed to detect main ocean is: ', jiseed, jjseed
126         WRITE(numout,'(a,f8.0)' )'    distance (m) between namelist seed location and model seed location: ', zdistseed
127         WRITE(numout,*)
128      END IF
129     
130      WHERE (msk_opnsea(:,:) > 0) msk_opnsea(:,:) = 0 ! mask all closed seas
131      WHERE (msk_opnsea(:,:) < 0) msk_opnsea(:,:) = 1 ! restore mask value
132
133      !! 1.4 Define closed sea mask (all of them, ie defined in the namelist or not)
134      !! needed to remove the undefined closed seas at the end
135      msk_closea = (ssmask - msk_opnsea) * 99
136
137      !!----------------------------------------------------------------------
138      !! 2 : compute closed sea mask for global, rnf and emp cases
139      !!----------------------------------------------------------------------
140
141      !! 2.1 : initialisation of masks and lake indexes
142      jglo = 1 ; jrnf = 1 ; jemp = 1
143      !! mask used to group lake by net evap/precip distribution technics
144      msk_glo = msk_closea
145      msk_rnf = msk_closea
146      msk_emp = msk_closea
147
148      !! mask used to group multiple lake with the same river mouth (great lake for example)
149      msk_gloid = 0.0_wp
150      msk_rnfid = 0.0_wp
151      msk_empid = 0.0_wp
152
153      IF (lwp) WRITE(numout,*)'closed seas: '
154 
155      DO jcs = 1,nn_closea
156
157         !! 2.2 : print out and sanity check
158         !! flag changed in sanity checks
159         lskip = .FALSE.
160
161         !! define seed to extract the closed sea jcs
162         CALL dom_ngb(sn_lake(jcs)%rlonsrc, sn_lake(jcs)%rlatsrc, jiseed, jjseed, zdistseed, 'T')
163
164         !! how the excess of the closed seas is spread out:
165         IF (lwp) THEN
166            WRITE(numout,*)
167            WRITE(numout,'(a,a10,a,2f7.2,a)')'    processing lake ',TRIM(sn_lake(jcs)%cname)             &
168                 &        ,' ( lat/lon : ',sn_lake(jcs)%rlatsrc, sn_lake(jcs)%rlonsrc,' )'
169         END IF
170         cloc = sn_lake(jcs)%cloctrg
171         csch = sn_lake(jcs)%cschtrg
172 
173         !! set up indexes and mask  ! SELECT CASE ...
174         SELECT CASE (csch)
175         CASE('glo')
176            jsch = jglo
177            zmsksrc = msk_glo
178            zmsktrg = msk_gloid
179            IF (lwp) WRITE(numout,*)'        net evap/precip will be spread globally'
180         CASE('rnf') 
181            jsch = jrnf
182            zmsksrc = msk_rnf
183            zmsktrg = msk_rnfid
184            IF (lwp) WRITE(numout,*)'        net precip will be spread locally and net evap globally'
185         CASE('emp')
186            jsch = jemp 
187            zmsksrc = msk_emp
188            zmsktrg = msk_empid
189            IF (lwp) WRITE(numout,*)'        net precip will be spread locally'
190         CASE DEFAULT
191            CALL ctl_stop( 'STOP', 'domclo: unknown target type for lake (should be glo, emp, rnf)' )
192         END SELECT
193
194         !! sanity check on the msk value (closea mask should be 99 otherwise, lake already processed)
195         !! check if some lake are not connected
196         zchk = 0._wp
197         IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) zchk = msk_closea(mi0(jiseed),mj0(jjseed))
198         CALL mpp_max('domclo',zchk)
199         
200         IF (zchk(1) /= 99._wp) THEN
201            IF (lwp) THEN
202               WRITE(numout,*)
203               WRITE(numout,*)'W A R N I N G: lake ',TRIM(sn_lake(jcs)%cname),' already processed or seed on land, we skip it.'
204               WRITE(numout,*)'It is possible this lake is connected to another one already processed.'
205               WRITE(numout,*)'It is possible the seed is on land because the lake is not present in the bathymetry file'
206               WRITE(numout,*)'If concerned about it, check input file and namelist definition of this lake.'
207            END IF
208            lskip = .TRUE.
209         END IF
210
211         !! 2.3 : compute mask for global, rnf, emp case
212         IF (.NOT. lskip) THEN
213            !! fill close sea mask with counter value
214            CALL fill_pool( jiseed, jjseed, zmsksrc, REAL(jsch  ,8))
215            WHERE (zmsksrc == REAL(jsch,8))
216               msk_closea = 0._wp
217               zmsktrg    = sn_lake(jcs)%idtrg
218            END WHERE
219
220            !! compute location of river mouth and distance to river mouth
221            IF (cloc /= 'global') THEN
222               !! set a minimum value for radius of the river influence
223               zradtrg = 0._wp
224               IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) THEN
225                  zradtrg(1) = MAX( sn_lake(jcs)%radtrg, e1t( mi0(jiseed), mj0(jjseed) ), e2t( mi0(jiseed), mj0(jjseed) ) )
226               END IF
227               CALL mpp_max('domclo',zradtrg)
228             
229               !! compute seed location for print
230               CALL dom_ngb(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, jiseed, jjseed, zdistseed, 'T')
231               !! compute distance to river mouth
232               zdist(:,:) = dist(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, glamt, gphit)
233            END IF
234
235            !! define estuary
236            !! deal with local/coast cases
237            SELECT CASE (cloc)
238            CASE ('global')
239               WHERE (msk_opnsea(:,:) == 1) zmsktrg = sn_lake(jcs)%idtrg
240
241            CASE ('local')
242               !! compute mask
243               WHERE (zdist(:,:) < zradtrg(1) .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
244               !! print
245               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a)')'         river mouth area is defined by         points within ',zradtrg(1)         &
246                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                                  &
247                             &        ,' (closest point is i/j ',jiseed, jjseed,')'
248 
249            CASE ('coast')
250               !! define coastline mask
251               zmsk_coastline = 0._wp
252               DO jj=2,jpj-1
253                  DO ji=2,jpi-1
254                     IF ( ssmask(ji,jj) == 1._wp .AND. SUM(ssmask(ji-1:ji+1,jj-1:jj+1)) < 9 ) zmsk_coastline(ji,jj) = 1._wp
255                  END DO
256               END DO
257               CALL lbc_lnk('domclo', zmsk_coastline,'T',1._wp)
258   
259               !! compute mask
260               WHERE ( zdist(:,:) < zradtrg(1) .AND. zmsk_coastline(:,:) == 1 .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
261               !! print
262               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a)')'         river mouth area is defined by coastal points within ',zradtrg(1) &
263                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                                  &
264                             &        ,' (closest point is i/j ',jiseed, jjseed,')'
265
266            CASE DEFAULT
267               CALL ctl_stop( 'STOP', 'domclo: unknown target type for lake (should be global, coast or local)' )
268
269            END SELECT
270   
271            !! sanity check
272            zarea = glob_sum('domclo',zmsktrg * msk_opnsea)
273            IF (zarea == 0._wp) CALL ctl_stop( 'STOP', 'river mouth area is 0, tune lon/lat trg or radtrg for this lake') 
274
275            !! set up indexes and mask
276            SELECT CASE (csch)
277            CASE ('glo')
278               jglo = jglo + 1
279               msk_glo   = zmsksrc
280               msk_gloid = zmsktrg
281            CASE ('rnf')
282               jrnf = jrnf + 1
283               msk_rnf   = zmsksrc
284               msk_rnfid = zmsktrg
285            CASE ('emp')
286               jemp = jemp + 1
287               msk_emp   = zmsksrc
288               msk_empid = zmsktrg
289            END SELECT
290
291         END IF ! lskip
292
293      END DO ! nn_closea
294
295      !!----------------------------------------------------------------------
296      !! 3 : clean the masks of possible remaining undefined closed seas
297      !!----------------------------------------------------------------------
298
299      !! mask all the cells not defined as closed sea
300      WHERE ( msk_glo == 99 ) msk_glo = 0
301      WHERE ( msk_rnf == 99 ) msk_rnf = 0
302      WHERE ( msk_emp == 99 ) msk_emp = 0
303
304      !!  non defined closed sea
305      WHERE (msk_closea > 0) msk_closea = 1
306
307   END SUBROUTINE dom_clo
308
309   !!======================================================================
310
311END MODULE domclo
312
Note: See TracBrowser for help on using the repository browser.