[11201] | 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 |
---|
[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 | !!---------------------------------------------------------------------- |
---|
| 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 |
---|
[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] | 97 | 901 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] | 100 | 902 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 | |
---|
| 329 | END MODULE domclo |
---|
| 330 | |
---|