[11308] | 1 | MODULE domisf |
---|
| 2 | !!============================================================================== |
---|
| 3 | !! *** MODULE domisf *** |
---|
| 4 | !! Ocean domain : |
---|
| 5 | !!============================================================================== |
---|
| 6 | !! History : 4.1 ! 2019-07 (P. Mathiot) re-write of the geometry definition under ice shelf |
---|
| 7 | !!---------------------------------------------------------------------- |
---|
| 8 | !!---------------------------------------------------------------------- |
---|
| 9 | !! zgr_isf_nam : read ice shelf namelist; print and compatibility option check |
---|
| 10 | !! zgr_isf_zspace : check compatibility between bathymetry and ice shelf draft and adjust ice shelf draft if needed (z space) |
---|
| 11 | !! zgr_isf_kspace : adjust ice shelf draft and compute top level to fit 2 cell in the water column (k space) |
---|
| 12 | !! zgr_isf : compute top partial cell scale factor |
---|
| 13 | !! zgr_isf_e3uv_w : correct e3uw and e3vw in case of 2 cell in water column under an ice shelf |
---|
| 14 | !!--------------------------------------------------------------------- |
---|
| 15 | USE dom_oce |
---|
[13204] | 16 | USE domutl ! flood filling algorithm |
---|
[11308] | 17 | USE domngb ! find nearest neighbourg |
---|
| 18 | USE in_out_manager ! I/O manager |
---|
| 19 | USE lbclnk ! lbclnk and lbclnk_multi |
---|
| 20 | USE lib_mpp ! |
---|
| 21 | |
---|
| 22 | IMPLICIT NONE |
---|
| 23 | PRIVATE |
---|
| 24 | |
---|
| 25 | PUBLIC zgr_isf_nam !: read domisf namelist |
---|
| 26 | PUBLIC zgr_isf_zspace !: check compatibility ice shelf draft/bathymetry and dig ice shelf draft if needed |
---|
| 27 | PUBLIC zgr_isf_kspace !: compute misfdep and dig ice shelf draft if needed |
---|
| 28 | PUBLIC zps_isf !: compute zps scale factor |
---|
| 29 | PUBLIC zps_isf_e3uv_w !: compute e3uw and e3vw in case of 2 cell in water column |
---|
| 30 | PUBLIC zgr_isf_subgl |
---|
| 31 | |
---|
| 32 | INTEGER :: nn_kisfmax = 999 !: maximal number of level change allowed by ln_isfconnect option |
---|
| 33 | REAL(wp), PUBLIC :: rn_isfdep_min = 10.0_wp !: ice shelf minimal thickness |
---|
| 34 | REAL(wp), PUBLIC :: rn_glhw_min = 1.0e-3 !: threshold on hw to define grounding line water |
---|
| 35 | REAL(wp), PUBLIC :: rn_isfhw_min = 1.0e-3 !: threshold on hw to define isf draft into the cavity |
---|
| 36 | REAL(wp), PUBLIC :: rn_isfshallow = 0.0_wp !: threshold to define shallow ice shelf cavity |
---|
| 37 | REAL(wp), PUBLIC :: rn_zisfmax = 6000.0_wp !: maximun meter of ice we are allowed to dig to assure connectivity |
---|
| 38 | REAL(wp), PUBLIC :: rn_isfsubgllon, rn_isfsubgllat !: seed for the detection of the 3d open ocean |
---|
| 39 | LOGICAL , PUBLIC :: ln_isfcheminey= .FALSE. !: remove cheminey |
---|
| 40 | LOGICAL , PUBLIC :: ln_isfconnect = .FALSE. !: assure connectivity |
---|
| 41 | LOGICAL , PUBLIC :: ln_isfchannel = .FALSE. !: remove channel in water column (on z space not level space) |
---|
| 42 | LOGICAL , PUBLIC :: ln_isfsubgl = .FALSE. !: remove subglacial lakes |
---|
| 43 | |
---|
| 44 | CONTAINS |
---|
| 45 | |
---|
| 46 | SUBROUTINE zgr_isf_nam |
---|
| 47 | !!---------------------------------------------------------------------- |
---|
| 48 | !! *** ROUTINE zps_isf *** |
---|
| 49 | !! |
---|
| 50 | !! ** Purpose : Read namzgr_isf namelist and |
---|
| 51 | !! check compatibility with other option |
---|
| 52 | !! |
---|
| 53 | !!---------------------------------------------------------------------- |
---|
| 54 | !! |
---|
| 55 | INTEGER :: ios, ierr |
---|
| 56 | !!--------------------------------------------------------------------- |
---|
| 57 | NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, rn_isfdep_min, rn_glhw_min, rn_isfhw_min, & |
---|
| 58 | & ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat |
---|
| 59 | ! |
---|
| 60 | ! 0.0 read namelist |
---|
[14623] | 61 | ! REWIND( numnam_ref ) ! Namelist namzgr_isf in reference namelist : ice shelf geometry definition |
---|
[11308] | 62 | READ ( numnam_ref, namzgr_isf, IOSTAT = ios, ERR = 901) |
---|
[14623] | 63 | 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in reference namelist') |
---|
[11308] | 64 | |
---|
[14623] | 65 | ! REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : ice shelf geometry definition |
---|
[11308] | 66 | READ ( numnam_cfg, namzgr_isf, IOSTAT = ios, ERR = 902 ) |
---|
[14623] | 67 | 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in configuration namelist') |
---|
[11308] | 68 | IF(lwm) WRITE ( numond, namzgr_isf ) |
---|
| 69 | ! |
---|
| 70 | IF (lwp) THEN |
---|
| 71 | WRITE(numout,*) |
---|
| 72 | WRITE(numout,*) ' dom_isf : ' |
---|
| 73 | WRITE(numout,*) ' ~~~~~~~ ' |
---|
| 74 | WRITE(numout,*) ' nn_kisfmax = ',nn_kisfmax ! |
---|
| 75 | WRITE(numout,*) ' rn_zisfmax = ',rn_zisfmax ! |
---|
| 76 | WRITE(numout,*) ' rn_isfdep_min = ',rn_isfdep_min ! |
---|
| 77 | WRITE(numout,*) ' rn_isfhw_min = ',rn_isfhw_min ! |
---|
| 78 | WRITE(numout,*) ' rn_glhw_min = ',rn_glhw_min ! |
---|
| 79 | WRITE(numout,*) ' ln_isfcheminey = ',ln_isfcheminey ! |
---|
| 80 | WRITE(numout,*) ' ln_isfconnect = ',ln_isfconnect ! |
---|
| 81 | WRITE(numout,*) ' ln_isfchannel = ',ln_isfchannel ! |
---|
| 82 | END IF |
---|
| 83 | ! |
---|
| 84 | ! 0.1 compatibility option |
---|
| 85 | ierr = 0 |
---|
| 86 | IF ( ln_zco .OR. ln_sco ) ierr = ierr + 1 |
---|
| 87 | IF ( ierr > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' ) |
---|
| 88 | |
---|
| 89 | END SUBROUTINE zgr_isf_nam |
---|
| 90 | |
---|
| 91 | SUBROUTINE zgr_isf_zspace |
---|
| 92 | !!---------------------------------------------------------------------- |
---|
| 93 | !! *** ROUTINE dom_isf *** |
---|
| 94 | !! |
---|
| 95 | !! ** Purpose : compute the grounding line position |
---|
| 96 | !! |
---|
| 97 | !! ** Method : bathy set to 0 and isfdraft are modified (dig/close) to respect |
---|
| 98 | !! these criterions. |
---|
| 99 | !! digging and filling base on depth criterion only |
---|
| 100 | !! 1.0 = set iceshelf to the minimum depth allowed |
---|
| 101 | !! 1.1 = ground ice shelf if water column less than X m |
---|
| 102 | !! 1.2 = ensure a minimum thickness for iceshelf cavity in shallow water |
---|
| 103 | !! 1.3 = remove channels and single point 'bay' |
---|
| 104 | !! |
---|
| 105 | !! ** Action : - test compatibility between isfdraft and bathy |
---|
| 106 | !! - bathy and isfdraft are modified |
---|
| 107 | !!---------------------------------------------------------------------- |
---|
| 108 | !! |
---|
| 109 | INTEGER :: ji, jj ! loop indexes |
---|
| 110 | INTEGER :: imskjp1, imskjm1, imskip1, imskim1 ! local variable |
---|
| 111 | INTEGER, DIMENSION(jpi,jpj) :: imask ! isf mask |
---|
| 112 | ! |
---|
| 113 | REAL(wp) :: zisfdep_min ! minimal ice shelf draft allowed |
---|
| 114 | !!--------------------------------------------------------------------- |
---|
| 115 | ! |
---|
| 116 | ! 1.0 set iceshelf to the minimum depth allowed |
---|
| 117 | zisfdep_min = MAX(rn_isfdep_min,e3t_1d(1)) |
---|
| 118 | WHERE(risfdep(:,:) > 0.0_wp .AND. risfdep(:,:) < zisfdep_min) |
---|
| 119 | risfdep(:,:) = zisfdep_min |
---|
| 120 | END WHERE |
---|
| 121 | ! |
---|
| 122 | ! 1.1 ground ice shelf if water column less than rn_glhw_min m |
---|
| 123 | ! => set the grounding line position |
---|
[11568] | 124 | WHERE( bathy(:,:) - risfdep(:,:) < rn_glhw_min .AND. risfdep(:,:) > 0.0_wp ) |
---|
[11308] | 125 | bathy (:,:) = 0._wp ; risfdep(:,:) = 0._wp |
---|
| 126 | END WHERE |
---|
| 127 | ! |
---|
| 128 | ! 1.2 ensure a minimum thickness for iceshelf cavity |
---|
| 129 | ! => avoid to negative e3t if ssh + sum(e3t_0*tmask) < 0 |
---|
| 130 | DO jj = 1, jpj |
---|
| 131 | DO ji = 1, jpi |
---|
[11568] | 132 | IF ( bathy(ji,jj) - risfdep(ji,jj) < rn_isfhw_min .AND. risfdep(ji,jj) > 0.0_wp ) THEN |
---|
[11308] | 133 | risfdep(ji,jj) = bathy(ji,jj) - rn_isfhw_min |
---|
| 134 | ! sanity check on risfdep (if < zisfdep_min) |
---|
| 135 | ! => we ground it as it failed to respect condition 1.0 and 1.1 |
---|
| 136 | IF ( risfdep(ji,jj) < zisfdep_min ) THEN |
---|
| 137 | bathy(ji,jj)=0._wp ; risfdep(ji,jj)=0._wp |
---|
| 138 | END IF |
---|
| 139 | END IF |
---|
| 140 | END DO |
---|
| 141 | END DO |
---|
| 142 | ! |
---|
| 143 | ! 1.3 Remove channels and single point 'bay' using bathy mask. |
---|
| 144 | ! => channel could be created if connectivity is enforced later. |
---|
| 145 | IF (ln_isfchannel) THEN |
---|
| 146 | imask(:,:) = 0 |
---|
[11568] | 147 | WHERE ( bathy(:,:) > 0._wp ) |
---|
[11308] | 148 | imask(:,:) = 1 |
---|
| 149 | END WHERE |
---|
| 150 | DO jj = 2, jpjm1 |
---|
| 151 | DO ji = 2, jpim1 |
---|
[11568] | 152 | IF( imask(ji,jj) == 1 .AND. risfdep(ji,jj) > 0._wp) THEN |
---|
[11308] | 153 | ! define connexion |
---|
| 154 | imskip1=imask(ji,jj)*imask(ji+1,jj ) ! 1 = connexion |
---|
| 155 | imskjp1=imask(ji,jj)*imask(ji ,jj+1) ! 1 = connexion |
---|
| 156 | imskim1=imask(ji,jj)*imask(ji-1,jj ) ! 1 = connexion |
---|
| 157 | imskjm1=imask(ji,jj)*imask(ji ,jj-1) ! 1 = connexion |
---|
| 158 | ! zonal channel and single bay |
---|
| 159 | IF ((imskip1+imskim1>=1) .AND. (imskjp1+imskjm1==0)) THEN |
---|
| 160 | bathy(ji,jj)=0._wp ; risfdep(ji,jj)=0._wp |
---|
| 161 | END IF |
---|
| 162 | ! meridionnal channel and single bay |
---|
| 163 | IF ((imskjp1+imskjm1>=1) .AND. (imskip1+imskim1==0)) THEN |
---|
| 164 | bathy(ji,jj)=0._wp ; risfdep(ji,jj)=0._wp |
---|
| 165 | END IF |
---|
| 166 | END IF |
---|
| 167 | END DO |
---|
| 168 | END DO |
---|
| 169 | ! ensure halo correct |
---|
| 170 | CALL lbc_lnk_multi( 'domisf', risfdep, 'T', 1._wp, bathy , 'T', 1._wp ) |
---|
| 171 | END IF |
---|
| 172 | ! |
---|
| 173 | END SUBROUTINE zgr_isf_zspace |
---|
| 174 | |
---|
| 175 | SUBROUTINE zgr_isf_kspace |
---|
| 176 | !!---------------------------------------------------------------------- |
---|
| 177 | !! *** ROUTINE dom_isf *** |
---|
| 178 | !! |
---|
| 179 | !! ** Purpose : compute misfdep and update isf draft if needed |
---|
| 180 | !! |
---|
| 181 | !! ** Method : The water column has to contain at least 2 wet cells in the water column under an ice shelf |
---|
| 182 | !! isf draft is modified (dig/close) to respect this criterion. |
---|
| 183 | !! compute level |
---|
| 184 | !! 2.0 = compute level |
---|
| 185 | !! 2.1 = ensure misfdep is not below bathymetry after step 2.0 |
---|
| 186 | !! digging and filling base on level criterion only |
---|
| 187 | !! 3.0 = dig to fit the 2 water level criterion (closed pool possible after this step) |
---|
| 188 | !! 3.1 = dig enough to ensure connectivity of all the cell beneath an ice shelf |
---|
| 189 | !! (most of the close pool should be remove after this step) |
---|
| 190 | !! 3.2 = fill chimney |
---|
| 191 | !! |
---|
| 192 | !! ** Action : - compute misfdep |
---|
| 193 | !! - isf draft is modified if needed |
---|
| 194 | !!---------------------------------------------------------------------- |
---|
| 195 | INTEGER :: ji, jj, jk, jn |
---|
| 196 | INTEGER :: icompt |
---|
| 197 | INTEGER :: ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 |
---|
| 198 | INTEGER :: ibathy, ibathyij, ibathyim1, ibathyip1, ibathyjm1, ibathyjp1 |
---|
| 199 | INTEGER :: imskjp1 , imskjm1 , imskip1 , imskim1 |
---|
| 200 | INTEGER :: imskip1_r, imskim1_r, imskjp1_r, imskjm1_r |
---|
| 201 | INTEGER , DIMENSION(jpi,jpj) :: imbathy, imisfdep, imask |
---|
| 202 | ! |
---|
| 203 | REAL(wp) :: zdepth |
---|
| 204 | REAL(wp), DIMENSION(jpi,jpj) :: zrisfdep, zdummy ! 2D workspace (ISH) |
---|
| 205 | !!--------------------------------------------------------------------- |
---|
| 206 | ! |
---|
| 207 | ! 2 Compute misfdep for ocean points (i.e. first wet level) |
---|
| 208 | ! find the first ocean level such that the first level thickness |
---|
| 209 | ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where |
---|
| 210 | ! e3t_0 is the reference level thickness |
---|
| 211 | ! |
---|
| 212 | ! 2.0 compute misfdep (taking into account the minimal cell thickness allowed) |
---|
| 213 | ! set misfdep to 1 on open water and initialisation beneath ice shelf |
---|
| 214 | WHERE( risfdep(:,:) == 0._wp ) ; misfdep(:,:) = 1 ! open water or grounded : set misfdep to 1 |
---|
| 215 | ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level |
---|
| 216 | END WHERE |
---|
| 217 | ! |
---|
| 218 | DO jk = 2, jpkm1 |
---|
| 219 | zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) |
---|
| 220 | WHERE( risfdep(:,:) > 0._wp .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 |
---|
| 221 | END DO |
---|
| 222 | ! |
---|
| 223 | ! 2.1 fill isolated grid point in the bathymetry |
---|
| 224 | ! will be done again later on in zgr_bat_ctl, but need to be done here to adjust misfdep respectively |
---|
| 225 | icompt = 0 |
---|
| 226 | DO jj = 2, jpjm1 |
---|
| 227 | DO ji = 2, jpim1 |
---|
| 228 | ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), & |
---|
| 229 | & mbathy(ji,jj-1), mbathy(ji,jj+1) ) |
---|
| 230 | IF( ibtest < mbathy(ji,jj) ) THEN |
---|
| 231 | mbathy(ji,jj) = ibtest |
---|
| 232 | icompt = icompt + 1 |
---|
| 233 | END IF |
---|
| 234 | END DO |
---|
| 235 | END DO |
---|
| 236 | ! |
---|
| 237 | ! ensure halo correct |
---|
| 238 | zdummy(:,:) = FLOAT( mbathy(:,:) ) ; CALL lbc_lnk('domisf', zdummy, 'T', 1._wp ) ; mbathy(:,:) = INT( zdummy(:,:) ) |
---|
| 239 | ! |
---|
| 240 | IF( lk_mpp ) CALL mpp_sum('domisf', icompt ) |
---|
| 241 | IF( icompt == 0 ) THEN |
---|
| 242 | IF(lwp) WRITE(numout,*)' no isolated ocean grid points' |
---|
| 243 | ELSE |
---|
| 244 | IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' |
---|
| 245 | ENDIF |
---|
| 246 | ! |
---|
| 247 | ! 2.2 be sure misfdep not below mbathy |
---|
| 248 | ! warning because of condition 4 we could have 'wet cell with misfdep below mbathy |
---|
| 249 | ! risfdep of these cells will be fix later on (see 3) |
---|
| 250 | WHERE( misfdep > mbathy ) misfdep(:,:) = MAX( 1, mbathy(:,:) ) |
---|
| 251 | ! |
---|
| 252 | ! 3.0 Assure 2 wet cells in the water column at T point and along the edge. |
---|
| 253 | ! find the deepest isfdep level that fit the 2 wet cell on the water column |
---|
| 254 | ! on all the sides (still need 4 pass) |
---|
| 255 | ! It need 4 iterations: if un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1 |
---|
| 256 | DO jn = 1, 4 |
---|
| 257 | imisfdep = misfdep |
---|
| 258 | DO jj = 2, jpjm1 |
---|
| 259 | DO ji = 2, jpim1 |
---|
| 260 | ! ISF cell only |
---|
| 261 | IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN |
---|
| 262 | ibathyij = mbathy(ji ,jj) |
---|
| 263 | ! |
---|
| 264 | ! set ground edge value to jpk to skip it later on |
---|
| 265 | ibathyip1 = mbathy(ji+1,jj) ; IF ( ibathyip1 < misfdep(ji,jj) ) ibathyip1 = jpk ! no wet cell in common on this edge |
---|
| 266 | ibathyim1 = mbathy(ji-1,jj) ; IF ( ibathyim1 < misfdep(ji,jj) ) ibathyim1 = jpk ! no wet cell in common on this edge |
---|
| 267 | ibathyjp1 = mbathy(ji,jj+1) ; IF ( ibathyjp1 < misfdep(ji,jj) ) ibathyjp1 = jpk ! no wet cell in common on this edge |
---|
| 268 | ibathyjm1 = mbathy(ji,jj-1) ; IF ( ibathyjm1 < misfdep(ji,jj) ) ibathyjm1 = jpk ! no wet cell in common on this edge |
---|
| 269 | ! |
---|
| 270 | ! find shallowest bathy level among the current cell and the neigbourging cells |
---|
| 271 | ibathy = MIN(ibathyij,ibathyip1,ibathyim1,ibathyjp1,ibathyjm1) |
---|
| 272 | ! |
---|
| 273 | ! update misfdep and risfdep if needed |
---|
| 274 | ! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the water column |
---|
| 275 | jk = MIN(misfdep(ji,jj),ibathy-1) |
---|
| 276 | IF ( jk < misfdep(ji,jj) ) THEN |
---|
| 277 | imisfdep(ji,jj) = jk |
---|
| 278 | risfdep(ji,jj) = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) |
---|
| 279 | END IF |
---|
| 280 | ENDIF |
---|
| 281 | END DO |
---|
| 282 | END DO |
---|
| 283 | misfdep=imisfdep |
---|
| 284 | ! |
---|
| 285 | ! ensure halo correct before new pass |
---|
| 286 | zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk('domisf', zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) |
---|
| 287 | CALL lbc_lnk('domisf', risfdep, 'T', 1. ) |
---|
| 288 | END DO ! jn |
---|
| 289 | ! |
---|
| 290 | ! 3.1 condition block to inssure connectivity everywhere beneath an ice shelf |
---|
| 291 | IF (ln_isfconnect) THEN |
---|
| 292 | imask(:,:) = 1 |
---|
| 293 | imbathy = mbathy |
---|
| 294 | imisfdep = misfdep |
---|
| 295 | zrisfdep = risfdep |
---|
| 296 | WHERE ( mbathy(:,:) == 0 ) |
---|
[14623] | 297 | imask(:,:) = 0 |
---|
[11308] | 298 | imbathy(:,:) = jpk |
---|
| 299 | END WHERE |
---|
| 300 | DO jj = 2, jpjm1 |
---|
| 301 | DO ji = 2, jpim1 |
---|
| 302 | IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN |
---|
| 303 | ! |
---|
[14623] | 304 | ! what it should be |
---|
[11308] | 305 | imskip1 = imask(ji,jj) * imask(ji+1,jj ) ! 1 = should be connected |
---|
| 306 | imskim1 = imask(ji,jj) * imask(ji-1,jj ) ! 1 = should be connected |
---|
| 307 | imskjp1 = imask(ji,jj) * imask(ji ,jj+1) ! 1 = should be connected |
---|
| 308 | imskjm1 = imask(ji,jj) * imask(ji ,jj-1) ! 1 = should be connected |
---|
| 309 | ! |
---|
[14623] | 310 | ! what it is |
---|
[11308] | 311 | imskip1_r=jpk ; imskim1_r=jpk; imskjp1_r=jpk; imskjm1_r=jpk |
---|
| 312 | IF (misfdep(ji,jj) > imbathy(ji+1,jj )) imskip1_r=1.0 ! 1 = no effective connection |
---|
| 313 | IF (misfdep(ji,jj) > imbathy(ji-1,jj )) imskim1_r=1.0 ! 1 = no effective connection |
---|
| 314 | IF (misfdep(ji,jj) > imbathy(ji ,jj+1)) imskjp1_r=1.0 ! 1 = no effective connection |
---|
| 315 | IF (misfdep(ji,jj) > imbathy(ji ,jj-1)) imskjm1_r=1.0 ! 1 = no effective connection |
---|
| 316 | ! |
---|
| 317 | ! defining level needed for connectivity |
---|
[14623] | 318 | ! imskip1 * imskip1_r == 1 means connections need to be enforce |
---|
[11308] | 319 | jk=MIN(imbathy(ji+1,jj ) * imskip1_r * imskip1, & |
---|
| 320 | & imbathy(ji-1,jj ) * imskim1_r * imskim1, & |
---|
| 321 | & imbathy(ji ,jj+1) * imskjp1_r * imskjp1, & |
---|
| 322 | & imbathy(ji ,jj-1) * imskjm1_r * imskjm1, & |
---|
[14623] | 323 | & jpk+1 ) ! add jpk in the MIN to avoid out of boundary later on |
---|
[11308] | 324 | ! |
---|
| 325 | ! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep |
---|
| 326 | imisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1) |
---|
| 327 | ! |
---|
| 328 | ! update isf draft if needed (need to be done here because next test is in meter and level) |
---|
| 329 | IF ( imisfdep(ji,jj) < misfdep(ji,jj) ) THEN |
---|
| 330 | jk = imisfdep(ji,jj) |
---|
| 331 | zrisfdep(ji,jj) = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) |
---|
| 332 | END IF |
---|
| 333 | |
---|
| 334 | ! sanity check |
---|
| 335 | ! check if we dig more than nn_kisfmax level or reach the surface |
---|
| 336 | ! check if we dig more than rn_zisfmax meter |
---|
| 337 | ! => if this is the case, undo what has been done before |
---|
| 338 | IF ( (misfdep(ji,jj)-imisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) & |
---|
| 339 | & .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj) )) ) THEN |
---|
| 340 | imisfdep(ji,jj)=misfdep(ji,jj) |
---|
| 341 | zrisfdep(ji,jj)=risfdep(ji,jj) |
---|
| 342 | END IF |
---|
| 343 | END IF |
---|
| 344 | END DO |
---|
| 345 | END DO |
---|
| 346 | misfdep=imisfdep |
---|
| 347 | risfdep=zrisfdep |
---|
| 348 | ! |
---|
| 349 | ! ensure halo correct |
---|
| 350 | zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk('domisf', zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) |
---|
| 351 | CALL lbc_lnk('domisf', risfdep, 'T', 1. ) |
---|
| 352 | END IF |
---|
| 353 | ! |
---|
| 354 | ! 3.2 fill hole in ice shelf (ie cell with no velocity point) |
---|
| 355 | ! => misfdep = MIN(misfdep at U, U-1, V, V-1) |
---|
| 356 | ! risfdep = gdepw(misfdep) (ie full cell) |
---|
| 357 | IF (ln_isfcheminey) THEN |
---|
| 358 | imisfdep = misfdep |
---|
| 359 | WHERE (mbathy == 0) imisfdep=jpk ! grounded |
---|
| 360 | DO jj = 2, jpjm1 |
---|
| 361 | DO ji = 2, jpim1 |
---|
| 362 | ibtest = imisfdep(ji ,jj ) |
---|
| 363 | ibtestim1 = imisfdep(ji-1,jj ) ; ibtestip1 = imisfdep(ji+1,jj ) |
---|
| 364 | ibtestjm1 = imisfdep(ji ,jj-1) ; ibtestjp1 = imisfdep(ji ,jj+1) |
---|
| 365 | ! |
---|
| 366 | ! correction in case isfdraft(ii,jj) deeper than bathy U-1/U+1/V-1/V+1 |
---|
| 367 | IF( imisfdep(ji,jj) > mbathy(ji-1,jj ) ) ibtestim1 = jpk ! mask at U-1 |
---|
| 368 | IF( imisfdep(ji,jj) > mbathy(ji+1,jj ) ) ibtestip1 = jpk ! mask at U+1 |
---|
| 369 | IF( imisfdep(ji,jj) > mbathy(ji ,jj-1) ) ibtestjm1 = jpk ! mask at V-1 |
---|
| 370 | IF( imisfdep(ji,jj) > mbathy(ji ,jj+1) ) ibtestjp1 = jpk ! mask at V+1 |
---|
| 371 | ! |
---|
| 372 | ! correction in case bathy(ii,jj) shallower than isfdraft U-1/U+1/V-1/V+1 |
---|
| 373 | IF( mbathy(ji,jj) < imisfdep(ji-1,jj ) ) ibtestim1 = jpk ! mask at U-1 |
---|
| 374 | IF( mbathy(ji,jj) < imisfdep(ji+1,jj ) ) ibtestip1 = jpk ! mask at U+1 |
---|
| 375 | IF( mbathy(ji,jj) < imisfdep(ji ,jj-1) ) ibtestjm1 = jpk ! mask at V-1 |
---|
| 376 | IF( mbathy(ji,jj) < imisfdep(ji ,jj+1) ) ibtestjp1 = jpk ! mask at V+1 |
---|
| 377 | ! |
---|
| 378 | ! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job) |
---|
| 379 | ! if misfdep is not changed, nothing is done |
---|
| 380 | ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1)) |
---|
| 381 | IF (misfdep(ji,jj) < ibtest) THEN |
---|
| 382 | misfdep(ji,jj) = ibtest |
---|
| 383 | risfdep(ji,jj) = gdepw_1d(ibtest) |
---|
| 384 | END IF |
---|
| 385 | ENDDO |
---|
| 386 | ENDDO |
---|
| 387 | ! |
---|
| 388 | ! if surround is fully masked, we mask it |
---|
| 389 | WHERE( misfdep(:,:) == jpk) ! ground case (1) |
---|
| 390 | mbathy(:,:) = 0 ; bathy(:,:) = 0.0_wp ; misfdep(:,:) = 1 ; risfdep(:,:) = 0.0_wp |
---|
| 391 | END WHERE |
---|
| 392 | ! |
---|
| 393 | ! ensure halo correct |
---|
| 394 | zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk('domisf', zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) ) |
---|
| 395 | zdummy(:,:) = FLOAT( mbathy (:,:) ); CALL lbc_lnk('domisf', zdummy, 'T', 1. ); mbathy (:,:) = INT( zdummy(:,:) ) |
---|
| 396 | CALL lbc_lnk('domisf', risfdep, 'T', 1. ) |
---|
| 397 | CALL lbc_lnk('domisf', bathy , 'T', 1. ) |
---|
| 398 | END IF |
---|
| 399 | ! |
---|
| 400 | END SUBROUTINE zgr_isf_kspace |
---|
| 401 | |
---|
| 402 | SUBROUTINE zgr_isf_subgl |
---|
| 403 | !!---------------------------------------------------------------------- |
---|
| 404 | !! *** ROUTINE zps_isf_subg *** |
---|
| 405 | !! |
---|
| 406 | !! ** Purpose : remove subglacial lakes |
---|
| 407 | !! |
---|
| 408 | !! ** Method : Use a flood filling algorithm to detect the open ocean |
---|
| 409 | !! and reverse the mask |
---|
| 410 | !! |
---|
| 411 | !! ** Action : - ckeck if seed on land |
---|
| 412 | !! - detect main ocean |
---|
| 413 | !! - mask subglacial lake (closed sea under an ice shelf) |
---|
| 414 | !! - mask properly mbathy, misfdep, bathy, risfdep |
---|
| 415 | !!---------------------------------------------------------------------- |
---|
| 416 | !! |
---|
| 417 | INTEGER :: jk ! loop index |
---|
| 418 | INTEGER :: iiseed, ijseed ! seed for flood filling algorithm |
---|
| 419 | INTEGER, DIMENSION(jpi,jpj) :: imask ! 2d mask |
---|
| 420 | ! |
---|
| 421 | REAL(wp) :: zdist ! distance between the seed and the nearest neighbourg |
---|
| 422 | REAL(wp), DIMENSION(1) :: zchk ! sanity check variable |
---|
| 423 | !!---------------------------------------------------------------------- |
---|
| 424 | ! |
---|
| 425 | ! check closed wet pool |
---|
| 426 | CALL dom_ngb(rn_isfsubgllon, rn_isfsubgllat, iiseed, ijseed, zdist, 'T') |
---|
| 427 | ! |
---|
| 428 | ! fill open ocean |
---|
| 429 | CALL fill_pool( iiseed, ijseed, 1, tmask, -1._wp ) |
---|
| 430 | ! at this point tmask:,:,:) can have 3 different values : |
---|
| 431 | ! 0 where there where already 0 |
---|
| 432 | ! -1 where the ocean points are connected |
---|
| 433 | ! 1 where ocean points in tmask are not connected |
---|
| 434 | IF (lwp) THEN |
---|
| 435 | WRITE(numout,*) |
---|
| 436 | WRITE(numout,*)'dom_isf_subg : removal of subglacial lakes ' |
---|
| 437 | WRITE(numout,*)'~~~~~~~~~~~~' |
---|
| 438 | WRITE(numout,*)'Number of disconected points : ', COUNT( (tmask(:,:,:) == 1) ) |
---|
| 439 | WRITE(numout,*)'lon/lat seed to detect main ocean is: ', rn_isfsubgllon, rn_isfsubgllat |
---|
| 440 | WRITE(numout,*)'i/j seed to detect main ocean is: ', iiseed, ijseed |
---|
| 441 | WRITE(numout,*) |
---|
| 442 | END IF |
---|
| 443 | ! |
---|
| 444 | DO jk = 1, jpk |
---|
| 445 | ! remove only subglacial lake (ie similar to close sea only below an ice shelf) |
---|
| 446 | WHERE (tmask(:,:,jk) > 0 .AND. misfdep(:,:) > 1) tmask(:,:,jk) = 0 |
---|
| 447 | WHERE (tmask(:,:,jk) < 0) tmask(:,:,jk) = 1 |
---|
| 448 | ! |
---|
| 449 | END DO |
---|
| 450 | ! |
---|
| 451 | ! surface mask |
---|
| 452 | imask(:,:) = MAXVAL( tmask(:,:,:), DIM=3 ) |
---|
| 453 | ! |
---|
| 454 | ! mask mbathy, misfdep, bathy and risfdep |
---|
| 455 | bathy(:,:) = bathy(:,:) * imask(:,:) |
---|
| 456 | mbathy(:,:) = mbathy(:,:) * imask(:,:) |
---|
| 457 | risfdep(:,:) = risfdep(:,:) * imask(:,:) |
---|
| 458 | misfdep(:,:) = misfdep(:,:) * imask(:,:) |
---|
| 459 | ! |
---|
| 460 | ! sanity check on seed location (land or not) |
---|
| 461 | zchk = 0._wp |
---|
| 462 | IF (mi0(iiseed) == mi1(iiseed) .AND. mj0(ijseed) == mj1(ijseed)) zchk = tmask(mi0(iiseed),mj0(ijseed),1) |
---|
| 463 | CALL mpp_max('domisf',zchk) |
---|
| 464 | IF (zchk(1) == 0._wp) CALL ctl_stop( 'STOP', 'open sea seed is on land, please update namelist (rn_lon_opnsea,rn_lat_opnsea)' ) |
---|
| 465 | ! |
---|
| 466 | END SUBROUTINE zgr_isf_subgl |
---|
| 467 | |
---|
| 468 | |
---|
| 469 | SUBROUTINE zps_isf |
---|
| 470 | !!---------------------------------------------------------------------- |
---|
| 471 | !! *** ROUTINE zps_isf *** |
---|
| 472 | !! |
---|
| 473 | !! ** Purpose : compute top partial cell |
---|
| 474 | !! |
---|
| 475 | !! ** Method : same method as for the bottom adapted for the top |
---|
| 476 | !! |
---|
| 477 | !!---------------------------------------------------------------------- |
---|
| 478 | !! |
---|
| 479 | INTEGER :: jj, ji ! loop indices |
---|
| 480 | INTEGER :: ik ! local top level |
---|
| 481 | INTEGER :: it ! counter index |
---|
| 482 | ! |
---|
| 483 | REAL(wp) :: zdiff |
---|
| 484 | !!---------------------------------------------------------------------- |
---|
| 485 | ! |
---|
| 486 | ! compute top partial cell |
---|
| 487 | DO jj = 1, jpj |
---|
| 488 | DO ji = 1, jpi |
---|
| 489 | ik = misfdep(ji,jj) |
---|
| 490 | IF( ik > 1 ) THEN ! ice shelf point only |
---|
| 491 | IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) |
---|
| 492 | gdepw_0(ji,jj,ik) = risfdep(ji,jj) |
---|
| 493 | !gm Bu check the gdepw_0 |
---|
| 494 | ! ... on ik |
---|
| 495 | gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & |
---|
| 496 | & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & |
---|
| 497 | & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) |
---|
| 498 | e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) |
---|
| 499 | e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) |
---|
| 500 | |
---|
| 501 | IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) |
---|
| 502 | e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) |
---|
| 503 | ENDIF |
---|
| 504 | ! ... on ik / ik-1 |
---|
| 505 | e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) |
---|
| 506 | e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) |
---|
| 507 | gdept_0(ji,jj,ik-1) = gdept_0(ji,jj,ik) - e3w_0(ji,jj,ik) |
---|
| 508 | gdepw_0(ji,jj,ik-1) = gdepw_0(ji,jj,ik) - e3t_0(ji,jj,ik-1) |
---|
[11568] | 509 | IF ( ik > 2 ) THEN |
---|
| 510 | e3w_0 (ji,jj,ik-1) = gdept_0(ji,jj,ik-1) - gdept_0(ji,jj,ik-2) |
---|
| 511 | ELSE |
---|
| 512 | e3w_0 (ji,jj,ik-1) = 2.0 * gdept_0(ji,jj,ik-1) |
---|
| 513 | END IF |
---|
[11308] | 514 | ENDIF |
---|
| 515 | END DO |
---|
| 516 | END DO |
---|
| 517 | ! |
---|
| 518 | it = 0 |
---|
| 519 | DO jj = 1, jpj |
---|
| 520 | DO ji = 1, jpi |
---|
| 521 | ik = misfdep(ji,jj) |
---|
| 522 | IF( ik > 1 ) THEN ! ice shelf point only |
---|
| 523 | e3tp (ji,jj) = e3t_0(ji,jj,ik ) |
---|
| 524 | e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) |
---|
| 525 | ! test |
---|
| 526 | zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) |
---|
| 527 | IF( zdiff <= 0. .AND. lwp ) THEN |
---|
| 528 | it = it + 1 |
---|
| 529 | WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj |
---|
| 530 | WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) |
---|
| 531 | WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff |
---|
| 532 | WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) |
---|
| 533 | ENDIF |
---|
| 534 | ENDIF |
---|
| 535 | END DO |
---|
| 536 | END DO |
---|
| 537 | ! |
---|
| 538 | END SUBROUTINE zps_isf |
---|
| 539 | |
---|
| 540 | SUBROUTINE zps_isf_e3uv_w |
---|
| 541 | !!---------------------------------------------------------------------- |
---|
| 542 | !! *** ROUTINE zps_isf *** |
---|
| 543 | !! |
---|
| 544 | !! ** Purpose : correct e3uw and e3vw for case with 2 wet cell in the water column under an ice shelf |
---|
| 545 | !! |
---|
| 546 | !!---------------------------------------------------------------------- |
---|
| 547 | !! |
---|
| 548 | INTEGER :: ji , jj ! loop indices |
---|
| 549 | INTEGER :: ikt, ikb ! local top/bottom level |
---|
| 550 | !!---------------------------------------------------------------------- |
---|
| 551 | ! |
---|
| 552 | ! define e3uw (adapted for 2 cells in the water column) |
---|
| 553 | DO jj = 2, jpjm1 |
---|
| 554 | DO ji = 2, jpim1 ! vector opt. |
---|
| 555 | ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) |
---|
| 556 | ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) |
---|
| 557 | IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & |
---|
| 558 | & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) |
---|
| 559 | ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) |
---|
| 560 | ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) |
---|
| 561 | IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & |
---|
| 562 | & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) |
---|
| 563 | END DO |
---|
| 564 | END DO |
---|
| 565 | ! |
---|
| 566 | END SUBROUTINE zps_isf_e3uv_w |
---|
| 567 | |
---|
| 568 | END MODULE domisf |
---|