source: utils/tools_ticket2457/DOMAINcfg/src/domisf.F90 @ 12871

Last change on this file since 12871 was 12871, checked in by mathiot, 9 months ago

ticket #2457: add fix suggested in ticket.

File size: 28.5 KB
Line 
1MODULE 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
44CONTAINS
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)
63901   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 )
67902   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
571END MODULE domisf
Note: See TracBrowser for help on using the repository browser.