[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 |
---|
| 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 |
---|
[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 ) |
---|
| 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 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 | |
---|
| 327 | END MODULE domclo |
---|
| 328 | |
---|