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

source: utils/tools/DOMAINcfg/src/domclo.F90

Last change on this file was 14623, checked in by ldebreu, 3 years ago

AGFdomcfg: 1) Update DOMAINcfg to be compliant with the removal of halo cells 2) Update most of the LBC ... subroutines to a recent NEMO 4 version #2638

File size: 15.5 KB
RevLine 
[11201]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
[13204]14   USE domutl          ! flood filling algorithm (fill_pool)
[11201]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
[11604]60         REAL(wp)       :: radtrg                      ! radius of closed sea river mouth (used if cschtrg is rnf or emp)
[11201]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
[12442]68      INTEGER :: jloc, jcoast       ! counter
[11201]69      INTEGER :: ios
70      INTEGER :: nn_closea          ! number of closed seas
71
[11604]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
[11201]78
[11604]79      CHARACTER(256) :: csch, cloc                       ! scheme name for water spreading (glo, rnf, emp)
80      TYPE(closea)  , DIMENSION(jp_closea)   :: sn_lake  ! lake properties
[11201]81
82      LOGICAL :: lskip     ! flag in case lake seed on land or already filled (...)
83
[13204]84      NAMELIST/namclo/ rn_lon_opnsea, rn_lat_opnsea, nn_closea, sn_lake
85
[11201]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      !!---------------------------------------------------------------------
[11604]94     
[14623]95   !   REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
[11201]96      READ  ( numnam_ref, namclo, IOSTAT = ios, ERR = 901 )
[14623]97901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namclo in reference namelist')
98   !   REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
[11201]99      READ  ( numnam_cfg, namclo, IOSTAT = ios, ERR = 902 )
[14623]100902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namclo in configuration namelist')
[11201]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
[11991]111      !! 1.2 fill connected cell to -1
[11201]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
[11604]121      !! print
[11201]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     
[11604]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
[11201]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
[11991]137      msk_csundef(:,:) = ( ssmask(:,:) - msk_opnsea(:,:) ) * 99._wp
[11201]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
[11991]146      msk_csglo(:,:) = msk_csundef(:,:)
147      msk_csrnf(:,:) = msk_csundef(:,:)
148      msk_csemp(:,:) = msk_csundef(:,:)
[11201]149
[11991]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
[11201]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
[11658]162
163         !! set up indexes and mask
[11201]164         SELECT CASE (csch)
165         CASE('glo')
166            jsch = jglo
[11991]167            zmsksrc(:,:) = msk_csglo(:,:)
168            zmsktrg(:,:) = msk_csgrpglo(:,:)
[11201]169         CASE('rnf') 
170            jsch = jrnf
[11991]171            zmsksrc(:,:) = msk_csrnf(:,:)
172            zmsktrg(:,:) = msk_csgrprnf(:,:)
[11201]173         CASE('emp')
174            jsch = jemp 
[11991]175            zmsksrc(:,:) = msk_csemp(:,:)
176            zmsktrg(:,:) = msk_csgrpemp(:,:)
[11201]177         CASE DEFAULT
[11658]178            CALL ctl_stop( 'STOP', 'domclo: ',TRIM(csch),' is an unknown target type for lake (should be glo, emp, rnf)' )
[11201]179         END SELECT
180
[11658]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
[11201]188         !! sanity check on the msk value (closea mask should be 99 otherwise, lake already processed)
[11991]189         !! check if seed on land or lake already processed
[11201]190         zchk = 0._wp
[11991]191         IF (mi0(jiseed) == mi1(jiseed) .AND. mj0(jjseed) == mj1(jjseed)) zchk = msk_csundef(mi0(jiseed),mj0(jjseed))
[11201]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
[11658]204 
[11201]205         !! 2.3 : compute mask for global, rnf, emp case
206         IF (.NOT. lskip) THEN
[11658]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
[11201]217            !! fill close sea mask with counter value
[11991]218            !! and update undefined closed sea mask
[11201]219            CALL fill_pool( jiseed, jjseed, zmsksrc, REAL(jsch  ,8))
[11604]220            WHERE (zmsksrc(:,:) == REAL(jsch,8))
[11991]221               msk_csundef = 0._wp
[11201]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
[11604]227
[11991]228               !! set value for radius of the river influence
[11201]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')
[11604]237
[11201]238               !! compute distance to river mouth
239               zdist(:,:) = dist(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, glamt, gphit)
[11604]240
[11201]241            END IF
242
243            !! define estuary
[11628]244            !! deal with global/local/coastal cases
[11201]245            SELECT CASE (cloc)
246            CASE ('global')
[11604]247               WHERE (msk_opnsea(:,:) == 1._wp) zmsktrg = sn_lake(jcs)%idtrg
[11201]248
249            CASE ('local')
250               !! compute mask
251               WHERE (zdist(:,:) < zradtrg(1) .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
[11604]252
[11201]253               !! print
[11658]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 )'
[11201]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
[11604]270
[11201]271               !! print
[11658]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 )'
[11201]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') 
[11658]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')
[11201]289
290            !! set up indexes and mask
291            SELECT CASE (csch)
292            CASE ('glo')
293               jglo = jglo + 1
[11991]294               msk_csglo(:,:)    = zmsksrc(:,:)
295               msk_csgrpglo(:,:) = zmsktrg(:,:)
[11658]296               IF (lwp) WRITE(numout,*)'        net evap/precip will be spread globally (glo)'
[11201]297            CASE ('rnf')
298               jrnf = jrnf + 1
[11991]299               msk_csrnf(:,:)    = zmsksrc(:,:)
300               msk_csgrprnf(:,:) = zmsktrg(:,:)
[11658]301               IF (lwp) WRITE(numout,*)'        net precip will be spread locally and net evap globally (rnf)'
[11201]302            CASE ('emp')
303               jemp = jemp + 1
[11991]304               msk_csemp(:,:)    = zmsksrc(:,:)
305               msk_csgrpemp(:,:) = zmsktrg(:,:)
[11658]306               IF (lwp) WRITE(numout,*)'        net evap/precip will be spread locally (emp)'
[11201]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
[11991]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
[11201]321
322      !!  non defined closed sea
[11991]323      WHERE ( msk_csundef(:,:) > 0._wp ) msk_csundef = 1._wp
[11201]324
325   END SUBROUTINE dom_clo
326
327   !!======================================================================
328
329END MODULE domclo
330
Note: See TracBrowser for help on using the repository browser.