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 utils/tools_ticket2457/DOMAINcfg/src – NEMO

source: utils/tools_ticket2457/DOMAINcfg/src/domclo.F90 @ 12871

Last change on this file since 12871 was 12442, checked in by mathiot, 4 years ago

changes required in tools after the merge of 2019 branche
File size: 15.5 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 (used if cschtrg is rnf or emp)
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 :: 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) :: zdist         ! distance to seed trg location
77      REAL(wp), DIMENSION(jpi,jpj) :: zmsksrc, zmsktrg, zmsk_coastline ! various mask
78
79      CHARACTER(256) :: csch, cloc                       ! scheme name for water spreading (glo, rnf, emp)
80      TYPE(closea)  , DIMENSION(jp_closea)   :: sn_lake  ! lake properties
81
82      LOGICAL :: lskip     ! flag in case lake seed on land or already filled (...)
83
84      !!----------------------------------------------------------------------
85      !! 0 : Read namelist for closed sea definition
86      !!----------------------------------------------------------------------
87      IF(lwp) WRITE(numout,*)
88      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
89      IF(lwp) WRITE(numout,*)'~~~~~~~'
90
91      NAMELIST/namclo/ rn_lon_opnsea, rn_lat_opnsea, nn_closea, sn_lake
92      !!---------------------------------------------------------------------
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 connected 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      !! print
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      !! 1.3 set to 0 everything >0 and revert mask
131      WHERE (msk_opnsea(:,:) > 0._wp) msk_opnsea(:,:) = 0._wp ! mask all closed seas
132      WHERE (msk_opnsea(:,:) < 0._wp) msk_opnsea(:,:) = 1._wp ! restore mask value
133
134      !! 1.4 Define closed sea mask (all of them, ie defined in the namelist or not)
135      !! needed to remove the undefined closed seas at the end
136      msk_csundef(:,:) = ( ssmask(:,:) - msk_opnsea(:,:) ) * 99._wp
137
138      !!----------------------------------------------------------------------
139      !! 2 : compute closed sea mask for global, rnf and emp cases
140      !!----------------------------------------------------------------------
141
142      !! 2.1 : initialisation of masks and lake indexes
143      jglo = 1 ; jrnf = 1 ; jemp = 1
144      !! mask used to group lake by net evap/precip distribution technics
145      msk_csglo(:,:) = msk_csundef(:,:)
146      msk_csrnf(:,:) = msk_csundef(:,:)
147      msk_csemp(:,:) = msk_csundef(:,:)
148
149      !! mask used to define group of lake sharing the same river mouth
150      msk_csgrpglo(:,:) = 0.0_wp
151      msk_csgrprnf(:,:) = 0.0_wp
152      msk_csgrpemp(:,:) = 0.0_wp
153
154      IF (lwp) WRITE(numout,*)'closed seas: '
155 
156      DO jcs = 1,nn_closea
157
158         !! how the excess of the closed seas is spread out:
159         cloc = sn_lake(jcs)%cloctrg
160         csch = sn_lake(jcs)%cschtrg
161
162         !! set up indexes and mask
163         SELECT CASE (csch)
164         CASE('glo')
165            jsch = jglo
166            zmsksrc(:,:) = msk_csglo(:,:)
167            zmsktrg(:,:) = msk_csgrpglo(:,:)
168         CASE('rnf') 
169            jsch = jrnf
170            zmsksrc(:,:) = msk_csrnf(:,:)
171            zmsktrg(:,:) = msk_csgrprnf(:,:)
172         CASE('emp')
173            jsch = jemp 
174            zmsksrc(:,:) = msk_csemp(:,:)
175            zmsktrg(:,:) = msk_csgrpemp(:,:)
176         CASE DEFAULT
177            CALL ctl_stop( 'STOP', 'domclo: ',TRIM(csch),' is an unknown target type for lake (should be glo, emp, rnf)' )
178         END SELECT
179
180         !! 2.2 : print out and sanity check
181         !! flag changed in sanity checks
182         lskip = .FALSE.
183
184         !! define seed to extract the closed sea jcs
185         CALL dom_ngb(sn_lake(jcs)%rlonsrc, sn_lake(jcs)%rlatsrc, jiseed, jjseed, zdistseed, 'T')
186
187         !! sanity check on the msk value (closea mask should be 99 otherwise, lake already processed)
188         !! check if seed on land or lake already processed
189         zchk = 0._wp
190         IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) zchk = msk_csundef(mi0(jiseed),mj0(jjseed))
191         CALL mpp_max('domclo',zchk)
192         
193         IF (zchk(1) /= 99._wp) THEN
194            IF (lwp) THEN
195               WRITE(numout,*)
196               WRITE(numout,*)'W A R N I N G: lake ',TRIM(sn_lake(jcs)%cname),' already processed or seed on land, we skip it.'
197               WRITE(numout,*)'It is possible this lake is connected to another one already processed.'
198               WRITE(numout,*)'It is possible the seed is on land because the lake is not present in the bathymetry file'
199               WRITE(numout,*)'If concerned about it, check input file and namelist definition of this lake.'
200            END IF
201            lskip = .TRUE.
202         END IF
203 
204         !! 2.3 : compute mask for global, rnf, emp case
205         IF (.NOT. lskip) THEN
206
207            IF (lwp) THEN
208               WRITE(numout,*)
209               WRITE(numout,'(a,a10,a,2f7.2,a,2i6,a,f7.2,a)')                               &
210                    &         '    processing lake ',TRIM(sn_lake(jcs)%cname)               &
211                    &        ,' ( lat/lon : ',sn_lake(jcs)%rlatsrc, sn_lake(jcs)%rlonsrc    &
212                    &        ,' , i/j     : ',jiseed,jjseed                                 &
213                    &        ,' , distance to seed : ',zdistseed,'m )'
214            END IF
215
216            !! fill close sea mask with counter value
217            !! and update undefined closed sea mask
218            CALL fill_pool( jiseed, jjseed, zmsksrc, REAL(jsch  ,8))
219            WHERE (zmsksrc(:,:) == REAL(jsch,8))
220               msk_csundef = 0._wp
221               zmsktrg    = sn_lake(jcs)%idtrg
222            END WHERE
223
224            !! compute location of river mouth and distance to river mouth
225            IF (cloc /= 'global') THEN
226
227               !! set value for radius of the river influence
228               zradtrg = 0._wp
229               IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) THEN
230                  zradtrg(1) = MAX( sn_lake(jcs)%radtrg, e1t( mi0(jiseed), mj0(jjseed) ), e2t( mi0(jiseed), mj0(jjseed) ) )
231               END IF
232               CALL mpp_max('domclo',zradtrg)
233             
234               !! compute seed location for print
235               CALL dom_ngb(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, jiseed, jjseed, zdistseed, 'T')
236
237               !! compute distance to river mouth
238               zdist(:,:) = dist(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, glamt, gphit)
239
240            END IF
241
242            !! define estuary
243            !! deal with global/local/coastal cases
244            SELECT CASE (cloc)
245            CASE ('global')
246               WHERE (msk_opnsea(:,:) == 1._wp) zmsktrg = sn_lake(jcs)%idtrg
247
248            CASE ('local')
249               !! compute mask
250               WHERE (zdist(:,:) < zradtrg(1) .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
251
252               !! print
253               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a,f7.0,a)')                                            &
254                             &         '         river mouth area is defined by         points within ',zradtrg(1) &
255                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                &
256                             &        ,' (closest wet point is i/j ',jiseed, jjseed,' at ',zdistseed,' m )'
257            CASE ('coast')
258               !! define coastline mask
259               zmsk_coastline = 0._wp
260               DO jj=2,jpj-1
261                  DO ji=2,jpi-1
262                     IF ( ssmask(ji,jj) == 1._wp .AND. SUM(ssmask(ji-1:ji+1,jj-1:jj+1)) < 9 ) zmsk_coastline(ji,jj) = 1._wp
263                  END DO
264               END DO
265               CALL lbc_lnk('domclo', zmsk_coastline,'T',1._wp)
266   
267               !! compute mask
268               WHERE ( zdist(:,:) < zradtrg(1) .AND. zmsk_coastline(:,:) == 1 .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
269
270               !! print
271               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a,f7.0,a)')                                            &
272                             &         '         river mouth area is defined by coastal points within ',zradtrg(1) &
273                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                &
274                             &        ,' (closest wet point is i/j ',jiseed, jjseed,' at ',zdistseed,' m )'
275
276            CASE DEFAULT
277               CALL ctl_stop( 'STOP', 'domclo: unknown target type for lake (should be global, coast or local)' )
278
279            END SELECT
280   
281            !! sanity check
282            zarea = glob_sum('domclo',zmsktrg * msk_opnsea)
283            IF (zarea == 0._wp) CALL ctl_stop( 'STOP', 'river mouth area is 0, tune lon/lat trg or radtrg for this lake') 
284            !
285            zarea = glob_sum('domclo',zmsksrc * msk_opnsea)
286            IF (zarea > 0._wp) CALL ctl_stop( 'STOP', 'closed seas and open ocean have common points, '                           &
287               &                                    , 'tune lon/lat src or check if your lake is really closed on the model grid')
288
289            !! set up indexes and mask
290            SELECT CASE (csch)
291            CASE ('glo')
292               jglo = jglo + 1
293               msk_csglo(:,:)    = zmsksrc(:,:)
294               msk_csgrpglo(:,:) = zmsktrg(:,:)
295               IF (lwp) WRITE(numout,*)'        net evap/precip will be spread globally (glo)'
296            CASE ('rnf')
297               jrnf = jrnf + 1
298               msk_csrnf(:,:)    = zmsksrc(:,:)
299               msk_csgrprnf(:,:) = zmsktrg(:,:)
300               IF (lwp) WRITE(numout,*)'        net precip will be spread locally and net evap globally (rnf)'
301            CASE ('emp')
302               jemp = jemp + 1
303               msk_csemp(:,:)    = zmsksrc(:,:)
304               msk_csgrpemp(:,:) = zmsktrg(:,:)
305               IF (lwp) WRITE(numout,*)'        net evap/precip will be spread locally (emp)'
306            END SELECT
307
308         END IF ! lskip
309
310      END DO ! nn_closea
311
312      !!----------------------------------------------------------------------
313      !! 3 : clean the masks of possible remaining undefined closed seas
314      !!----------------------------------------------------------------------
315
316      !! mask all the cells not defined as closed sea
317      WHERE ( msk_csglo(:,:) == 99._wp )  msk_csglo = 0._wp
318      WHERE ( msk_csrnf(:,:) == 99._wp )  msk_csrnf = 0._wp
319      WHERE ( msk_csemp(:,:) == 99._wp )  msk_csemp = 0._wp
320
321      !!  non defined closed sea
322      WHERE ( msk_csundef(:,:) > 0._wp ) msk_csundef = 1._wp
323
324   END SUBROUTINE dom_clo
325
326   !!======================================================================
327
328END MODULE domclo
329
Note: See TracBrowser for help on using the repository browser.