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

source: utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domclo.F90 @ 13024

Last change on this file since 13024 was 13024, checked in by rblod, 4 years ago

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