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, 4 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
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
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
[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
68      INTEGER :: jglo, jloc, jcoast ! counter
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
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      !!---------------------------------------------------------------------
[11604]93     
[11201]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
[11604]120      !! print
[11201]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     
[11604]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
[11201]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
[11604]136      msk_closea(:,:) = ( ssmask(:,:) - msk_opnsea(:,:) ) * 99._wp
[11201]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
[11604]145      msk_glo(:,:) = msk_closea(:,:)
146      msk_rnf(:,:) = msk_closea(:,:)
147      msk_emp(:,:) = msk_closea(:,:)
[11201]148
149      !! mask used to group multiple lake with the same river mouth (great lake for example)
[11604]150      msk_gloid(:,:) = 0.0_wp
151      msk_rnfid(:,:) = 0.0_wp
152      msk_empid(:,:) = 0.0_wp
[11201]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
[11658]161
162         !! set up indexes and mask
[11201]163         SELECT CASE (csch)
164         CASE('glo')
165            jsch = jglo
[11604]166            zmsksrc(:,:) = msk_glo(:,:)
167            zmsktrg(:,:) = msk_gloid(:,:)
[11201]168         CASE('rnf') 
169            jsch = jrnf
[11604]170            zmsksrc(:,:) = msk_rnf(:,:)
171            zmsktrg(:,:) = msk_rnfid(:,:)
[11201]172         CASE('emp')
173            jsch = jemp 
[11604]174            zmsksrc(:,:) = msk_emp(:,:)
175            zmsktrg(:,:) = msk_empid(:,:)
[11201]176         CASE DEFAULT
[11658]177            CALL ctl_stop( 'STOP', 'domclo: ',TRIM(csch),' is an unknown target type for lake (should be glo, emp, rnf)' )
[11201]178         END SELECT
179
[11658]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
[11201]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
[11658]203 
[11201]204         !! 2.3 : compute mask for global, rnf, emp case
205         IF (.NOT. lskip) THEN
[11658]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
[11201]216            !! fill close sea mask with counter value
217            CALL fill_pool( jiseed, jjseed, zmsksrc, REAL(jsch  ,8))
[11604]218            WHERE (zmsksrc(:,:) == REAL(jsch,8))
[11201]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
[11604]225
[11201]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')
[11604]235
[11201]236               !! compute distance to river mouth
237               zdist(:,:) = dist(sn_lake(jcs)%rlontrg, sn_lake(jcs)%rlattrg, glamt, gphit)
[11604]238
[11201]239            END IF
240
241            !! define estuary
[11628]242            !! deal with global/local/coastal cases
[11201]243            SELECT CASE (cloc)
244            CASE ('global')
[11604]245               WHERE (msk_opnsea(:,:) == 1._wp) zmsktrg = sn_lake(jcs)%idtrg
[11201]246
247            CASE ('local')
248               !! compute mask
249               WHERE (zdist(:,:) < zradtrg(1) .AND. msk_opnsea(:,:) == 1 ) zmsktrg = sn_lake(jcs)%idtrg
[11604]250
[11201]251               !! print
[11658]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 )'
[11201]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
[11604]268
[11201]269               !! print
[11658]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 )'
[11201]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') 
[11658]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')
[11201]287
288            !! set up indexes and mask
289            SELECT CASE (csch)
290            CASE ('glo')
291               jglo = jglo + 1
[11604]292               msk_glo(:,:)   = zmsksrc(:,:)
293               msk_gloid(:,:) = zmsktrg(:,:)
[11658]294               IF (lwp) WRITE(numout,*)'        net evap/precip will be spread globally (glo)'
[11201]295            CASE ('rnf')
296               jrnf = jrnf + 1
[11604]297               msk_rnf(:,:)   = zmsksrc(:,:)
298               msk_rnfid(:,:) = zmsktrg(:,:)
[11658]299               IF (lwp) WRITE(numout,*)'        net precip will be spread locally and net evap globally (rnf)'
[11201]300            CASE ('emp')
301               jemp = jemp + 1
[11604]302               msk_emp(:,:)   = zmsksrc(:,:)
303               msk_empid(:,:) = zmsktrg(:,:)
[11658]304               IF (lwp) WRITE(numout,*)'        net evap/precip will be spread locally (emp)'
[11201]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
[11604]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
[11201]319
320      !!  non defined closed sea
[11604]321      WHERE ( msk_closea(:,:) > 0._wp ) msk_closea = 1._wp
[11201]322
323   END SUBROUTINE dom_clo
324
325   !!======================================================================
326
327END MODULE domclo
328
Note: See TracBrowser for help on using the repository browser.