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