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 |
---|