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 @ 11628

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

ENHANCE-03_domcfg: add cpp_ENHANCE-03_domcfg.fcm and remove make_closea_masks.py (as now managed in the F90 code and remove compilation file (ticket #2143)
File size: 14.6 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         !! 2.2 : print out and sanity check
159         !! flag changed in sanity checks
160         lskip = .FALSE.
161
162         !! define seed to extract the closed sea jcs
163         CALL dom_ngb(sn_lake(jcs)%rlonsrc, sn_lake(jcs)%rlatsrc, jiseed, jjseed, zdistseed, 'T')
164
165         !! how the excess of the closed seas is spread out:
166         IF (lwp) THEN
167            WRITE(numout,*)
168            WRITE(numout,'(a,a10,a,2f7.2,a)')'    processing lake ',TRIM(sn_lake(jcs)%cname)             &
169                 &        ,' ( lat/lon : ',sn_lake(jcs)%rlatsrc, sn_lake(jcs)%rlonsrc,' )'
170         END IF
171         cloc = sn_lake(jcs)%cloctrg
172         csch = sn_lake(jcs)%cschtrg
173 
174         !! set up indexes and mask  ! SELECT CASE ...
175         SELECT CASE (csch)
176         CASE('glo')
177            jsch = jglo
178            zmsksrc(:,:) = msk_glo(:,:)
179            zmsktrg(:,:) = msk_gloid(:,:)
180            IF (lwp) WRITE(numout,*)'        net evap/precip will be spread globally'
181         CASE('rnf') 
182            jsch = jrnf
183            zmsksrc(:,:) = msk_rnf(:,:)
184            zmsktrg(:,:) = msk_rnfid(:,:)
185            IF (lwp) WRITE(numout,*)'        net precip will be spread locally and net evap globally'
186         CASE('emp')
187            jsch = jemp 
188            zmsksrc(:,:) = msk_emp(:,:)
189            zmsktrg(:,:) = msk_empid(:,:)
190            IF (lwp) WRITE(numout,*)'        net precip will be spread locally'
191         CASE DEFAULT
192            CALL ctl_stop( 'STOP', 'domclo: unknown target type for lake (should be glo, emp, rnf)' )
193         END SELECT
194
195         !! sanity check on the msk value (closea mask should be 99 otherwise, lake already processed)
196         !! check if some lake are not connected
197         zchk = 0._wp
198         IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) zchk = msk_closea(mi0(jiseed),mj0(jjseed))
199         CALL mpp_max('domclo',zchk)
200         
201         IF (zchk(1) /= 99._wp) THEN
202            IF (lwp) THEN
203               WRITE(numout,*)
204               WRITE(numout,*)'W A R N I N G: lake ',TRIM(sn_lake(jcs)%cname),' already processed or seed on land, we skip it.'
205               WRITE(numout,*)'It is possible this lake is connected to another one already processed.'
206               WRITE(numout,*)'It is possible the seed is on land because the lake is not present in the bathymetry file'
207               WRITE(numout,*)'If concerned about it, check input file and namelist definition of this lake.'
208            END IF
209            lskip = .TRUE.
210         END IF
211
212         !! 2.3 : compute mask for global, rnf, emp case
213         IF (.NOT. lskip) THEN
214            !! fill close sea mask with counter value
215            CALL fill_pool( jiseed, jjseed, zmsksrc, REAL(jsch  ,8))
216            WHERE (zmsksrc(:,:) == REAL(jsch,8))
217               msk_closea = 0._wp
218               zmsktrg    = sn_lake(jcs)%idtrg
219            END WHERE
220
221            !! compute location of river mouth and distance to river mouth
222            IF (cloc /= 'global') THEN
223
224               !! set a minimum value for radius of the river influence
225               zradtrg = 0._wp
226               IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) THEN
227                  zradtrg(1) = MAX( sn_lake(jcs)%radtrg, e1t( mi0(jiseed), mj0(jjseed) ), e2t( mi0(jiseed), mj0(jjseed) ) )
228               END IF
229               CALL mpp_max('domclo',zradtrg)
230             
231               !! compute seed location for print
232               CALL dom_ngb(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, jiseed, jjseed, zdistseed, 'T')
233
234               !! compute distance to river mouth
235               zdist(:,:) = dist(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, glamt, gphit)
236
237            END IF
238
239            !! define estuary
240            !! deal with global/local/coastal cases
241            SELECT CASE (cloc)
242            CASE ('global')
243               WHERE (msk_opnsea(:,:) == 1._wp) zmsktrg = sn_lake(jcs)%idtrg
244
245            CASE ('local')
246               !! compute mask
247               WHERE (zdist(:,:) < zradtrg(1) .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
248
249               !! print
250               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a)')'         river mouth area is defined by         points within ',zradtrg(1)         &
251                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                                  &
252                             &        ,' (closest point is i/j ',jiseed, jjseed,')'
253            CASE ('coast')
254               !! define coastline mask
255               zmsk_coastline = 0._wp
256               DO jj=2,jpj-1
257                  DO ji=2,jpi-1
258                     IF ( ssmask(ji,jj) == 1._wp .AND. SUM(ssmask(ji-1:ji+1,jj-1:jj+1)) < 9 ) zmsk_coastline(ji,jj) = 1._wp
259                  END DO
260               END DO
261               CALL lbc_lnk('domclo', zmsk_coastline,'T',1._wp)
262   
263               !! compute mask
264               WHERE ( zdist(:,:) < zradtrg(1) .AND. zmsk_coastline(:,:) == 1 .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
265
266               !! print
267               IF (lwp) WRITE(numout,'(a,f7.0,a,2f7.2,a,2i7,a)')'         river mouth area is defined by coastal points within ',zradtrg(1) &
268                             &        ,' m of lat/lon ', sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg                                  &
269                             &        ,' (closest point is i/j ',jiseed, jjseed,')'
270
271            CASE DEFAULT
272               CALL ctl_stop( 'STOP', 'domclo: unknown target type for lake (should be global, coast or local)' )
273
274            END SELECT
275   
276            !! sanity check
277            zarea = glob_sum('domclo',zmsktrg * msk_opnsea)
278            IF (zarea == 0._wp) CALL ctl_stop( 'STOP', 'river mouth area is 0, tune lon/lat trg or radtrg for this lake') 
279
280            !! set up indexes and mask
281            SELECT CASE (csch)
282            CASE ('glo')
283               jglo = jglo + 1
284               msk_glo(:,:)   = zmsksrc(:,:)
285               msk_gloid(:,:) = zmsktrg(:,:)
286            CASE ('rnf')
287               jrnf = jrnf + 1
288               msk_rnf(:,:)   = zmsksrc(:,:)
289               msk_rnfid(:,:) = zmsktrg(:,:)
290            CASE ('emp')
291               jemp = jemp + 1
292               msk_emp(:,:)   = zmsksrc(:,:)
293               msk_empid(:,:) = zmsktrg(:,:)
294            END SELECT
295
296         END IF ! lskip
297
298      END DO ! nn_closea
299
300      !!----------------------------------------------------------------------
301      !! 3 : clean the masks of possible remaining undefined closed seas
302      !!----------------------------------------------------------------------
303
304      !! mask all the cells not defined as closed sea
305      WHERE ( msk_glo(:,:) == 99._wp )  msk_glo = 0._wp
306      WHERE ( msk_rnf(:,:) == 99._wp )  msk_rnf = 0._wp
307      WHERE ( msk_emp(:,:) == 99._wp )  msk_emp = 0._wp
308
309      !!  non defined closed sea
310      WHERE ( msk_closea(:,:) > 0._wp ) msk_closea = 1._wp
311
312   END SUBROUTINE dom_clo
313
314   !!======================================================================
315
316END MODULE domclo
317
Note: See TracBrowser for help on using the repository browser.