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