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-03_domcfg/src – NEMO

source: NEMO/branches/2019/ENHANCE-03_domcfg/src/domclo.F90 @ 11658

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

ENHANCE-03_domcfg: forced local domain to read the overlap region + re-arange the ocean output print (ticket #2143)
File size: 15.4 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 :: 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) :: 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 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      !! 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_closea(:,:) = ( 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_glo(:,:) = msk_closea(:,:)
146      msk_rnf(:,:) = msk_closea(:,:)
147      msk_emp(:,:) = msk_closea(:,:)
148
149      !! mask used to group multiple lake with the same river mouth (great lake for example)
150      msk_gloid(:,:) = 0.0_wp
151      msk_rnfid(:,:) = 0.0_wp
152      msk_empid(:,:) = 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_glo(:,:)
167            zmsktrg(:,:) = msk_gloid(:,:)
168         CASE('rnf') 
169            jsch = jrnf
170            zmsksrc(:,:) = msk_rnf(:,:)
171            zmsktrg(:,:) = msk_rnfid(:,:)
172         CASE('emp')
173            jsch = jemp 
174            zmsksrc(:,:) = msk_emp(:,:)
175            zmsktrg(:,:) = msk_empid(:,:)
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 some lake are not connected
189         zchk = 0._wp
190         IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) zchk = msk_closea(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            CALL fill_pool( jiseed, jjseed, zmsksrc, REAL(jsch  ,8))
218            WHERE (zmsksrc(:,:) == REAL(jsch,8))
219               msk_closea = 0._wp
220               zmsktrg    = sn_lake(jcs)%idtrg
221            END WHERE
222
223            !! compute location of river mouth and distance to river mouth
224            IF (cloc /= 'global') THEN
225
226               !! set a minimum value for radius of the river influence
227               zradtrg = 0._wp
228               IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) THEN
229                  zradtrg(1) = MAX( sn_lake(jcs)%radtrg, e1t( mi0(jiseed), mj0(jjseed) ), e2t( mi0(jiseed), mj0(jjseed) ) )
230               END IF
231               CALL mpp_max('domclo',zradtrg)
232             
233               !! compute seed location for print
234               CALL dom_ngb(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, jiseed, jjseed, zdistseed, 'T')
235
236               !! compute distance to river mouth
237               zdist(:,:) = dist(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, glamt, gphit)
238
239            END IF
240
241            !! define estuary
242            !! deal with global/local/coastal cases
243            SELECT CASE (cloc)
244            CASE ('global')
245               WHERE (msk_opnsea(:,:) == 1._wp) zmsktrg = sn_lake(jcs)%idtrg
246
247            CASE ('local')
248               !! compute mask
249               WHERE (zdist(:,:) < zradtrg(1) .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
250
251               !! print
252               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a,f7.0,a)')                                            &
253                             &         '         river mouth area is defined by         points within ',zradtrg(1) &
254                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                &
255                             &        ,' (closest wet point is i/j ',jiseed, jjseed,' at ',zdistseed,' m )'
256            CASE ('coast')
257               !! define coastline mask
258               zmsk_coastline = 0._wp
259               DO jj=2,jpj-1
260                  DO ji=2,jpi-1
261                     IF ( ssmask(ji,jj) == 1._wp .AND. SUM(ssmask(ji-1:ji+1,jj-1:jj+1)) < 9 ) zmsk_coastline(ji,jj) = 1._wp
262                  END DO
263               END DO
264               CALL lbc_lnk('domclo', zmsk_coastline,'T',1._wp)
265   
266               !! compute mask
267               WHERE ( zdist(:,:) < zradtrg(1) .AND. zmsk_coastline(:,:) == 1 .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
268
269               !! print
270               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a,f7.0,a)')                                            &
271                             &         '         river mouth area is defined by coastal points within ',zradtrg(1) &
272                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                &
273                             &        ,' (closest wet point is i/j ',jiseed, jjseed,' at ',zdistseed,' m )'
274
275            CASE DEFAULT
276               CALL ctl_stop( 'STOP', 'domclo: unknown target type for lake (should be global, coast or local)' )
277
278            END SELECT
279   
280            !! sanity check
281            zarea = glob_sum('domclo',zmsktrg * msk_opnsea)
282            IF (zarea == 0._wp) CALL ctl_stop( 'STOP', 'river mouth area is 0, tune lon/lat trg or radtrg for this lake') 
283            !
284            zarea = glob_sum('domclo',zmsksrc * msk_opnsea)
285            IF (zarea > 0._wp) CALL ctl_stop( 'STOP', 'closed seas and open ocean have common points, '                           &
286               &                                    , 'tune lon/lat src or check if your lake is really closed on the model grid')
287
288            !! set up indexes and mask
289            SELECT CASE (csch)
290            CASE ('glo')
291               jglo = jglo + 1
292               msk_glo(:,:)   = zmsksrc(:,:)
293               msk_gloid(:,:) = zmsktrg(:,:)
294               IF (lwp) WRITE(numout,*)'        net evap/precip will be spread globally (glo)'
295            CASE ('rnf')
296               jrnf = jrnf + 1
297               msk_rnf(:,:)   = zmsksrc(:,:)
298               msk_rnfid(:,:) = zmsktrg(:,:)
299               IF (lwp) WRITE(numout,*)'        net precip will be spread locally and net evap globally (rnf)'
300            CASE ('emp')
301               jemp = jemp + 1
302               msk_emp(:,:)   = zmsksrc(:,:)
303               msk_empid(:,:) = zmsktrg(:,:)
304               IF (lwp) WRITE(numout,*)'        net evap/precip will be spread locally (emp)'
305            END SELECT
306
307         END IF ! lskip
308
309      END DO ! nn_closea
310
311      !!----------------------------------------------------------------------
312      !! 3 : clean the masks of possible remaining undefined closed seas
313      !!----------------------------------------------------------------------
314
315      !! mask all the cells not defined as closed sea
316      WHERE ( msk_glo(:,:) == 99._wp )  msk_glo = 0._wp
317      WHERE ( msk_rnf(:,:) == 99._wp )  msk_rnf = 0._wp
318      WHERE ( msk_emp(:,:) == 99._wp )  msk_emp = 0._wp
319
320      !!  non defined closed sea
321      WHERE ( msk_closea(:,:) > 0._wp ) msk_closea = 1._wp
322
323   END SUBROUTINE dom_clo
324
325   !!======================================================================
326
327END MODULE domclo
328
Note: See TracBrowser for help on using the repository browser.