New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domisf.F90 in utils/tools/DOMAINcfg/src – NEMO

source: utils/tools/DOMAINcfg/src/domisf.F90 @ 13204

Last change on this file since 13204 was 13204, checked in by smasson, 4 years ago

tools: update with tools_dev_r12970_AGRIF_CMEMS

File size: 28.3 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 domutl            ! 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, 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
124      WHERE( bathy(:,:) - risfdep(:,:) < rn_glhw_min .AND. risfdep(:,:) > 0.0_wp ) 
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
132            IF ( bathy(ji,jj) - risfdep(ji,jj) < rn_isfhw_min .AND. risfdep(ji,jj) > 0.0_wp ) THEN
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
147         WHERE ( bathy(:,:) > 0._wp )
148            imask(:,:) = 1
149         END WHERE
150         DO jj = 2, jpjm1
151            DO ji = 2, jpim1
152               IF( imask(ji,jj) == 1 .AND. risfdep(ji,jj) > 0._wp) THEN
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 )
297            imask(:,:) = 0
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                  !
304                  ! what it should be
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                  !
310                  ! what it is
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
318                  ! imskip1 * imskip1_r == 1 means connections need to be enforce
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, &
323                     &   jpk+1 ) ! add jpk in the MIN to avoid out of boundary later on
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)
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
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
568END MODULE domisf
Note: See TracBrowser for help on using the repository browser.