source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

File size: 19.1 KB
Line 
1MODULE iscplrst
2   !!======================================================================
3   !!                       ***  MODULE  iscplrst  ***
4   !! Ocean forcing: update the restart file in case of ice sheet/ocean coupling
5   !!=====================================================================
6   !! History :  NEMO  ! 2015-01 P. Mathiot: original
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   iscpl_stp          : step management
11   !!   iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
12   !!----------------------------------------------------------------------
13   USE dom_oce         ! ocean space and time domain
14   USE domwri          ! ocean space and time domain
15   USE domvvl, ONLY : dom_vvl_interpol
16   USE phycst          ! physical constants
17   USE sbc_oce         ! surface boundary condition variables
18   USE oce             ! global tra/dyn variable
19   USE in_out_manager  ! I/O manager
20   USE iom             ! I/O module
21   USE lib_mpp         ! MPP library
22   USE lib_fortran     ! MPP library
23   USE lbclnk          ! communication
24   USE iscplini        ! ice sheet coupling: initialisation
25   USE iscplhsb        ! ice sheet coupling: conservation
26
27   IMPLICIT NONE
28   PRIVATE
29   
30   PUBLIC   iscpl_stp          ! step management
31   !!
32   !! * Substitutions 
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE iscpl_stp
42      !!----------------------------------------------------------------------
43      !!                   ***  ROUTINE iscpl_stp  ***
44      !!
45      !! ** Purpose : compute initialisation
46      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
47      !!              compute correction term if needed
48      !!
49      !!----------------------------------------------------------------------
50      INTEGER  ::   inum0
51      REAL(wp), DIMENSION(jpi,jpj) ::   zsmask_b
52      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztmask_b, zumask_b, zvmask_b
53      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t_b  , ze3u_b  , ze3v_b 
54      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepw_b
55      CHARACTER(20) :: cfile
56      !!----------------------------------------------------------------------
57
58
59
60      !! get restart variable
61      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b   ) ! need to extrapolate T/S
62      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b   ) ! need to correct barotropic velocity
63      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b   ) ! need to correct barotropic velocity
64      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b   ) ! need to correct barotropic velocity
65      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:) )  ! need to compute temperature correction
66      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b(:,:,:) )  ! need to correct barotropic velocity
67      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b(:,:,:) )  ! need to correct barotropic velocity
68      CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:) ) ! need to interpol vertical profile (vvl)
69
70      !! read namelist
71      CALL iscpl_init()
72
73      !!  ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
74      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
75
76      !! compute correction if conservation needed
77      IF ( ln_hsb ) THEN
78         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
79         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
80      END IF
81         
82      !! print mesh/mask
83      IF( nn_msh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
84
85      IF ( ln_hsb ) THEN
86         cfile='correction'
87         cfile = TRIM( cfile )
88         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
89         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
90         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
91         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
92         CALL iom_close ( inum0 )
93      END IF
94
95
96      !! next step is an euler time step
97      neuler = 0
98
99      !! set _b and _n variables equal
100      tsb (:,:,:,:) = tsn (:,:,:,:)
101      ub  (:,:,:)   = un  (:,:,:)
102      vb  (:,:,:)   = vn  (:,:,:)
103      sshb(:,:)     = sshn(:,:)
104
105      !! set _b and _n vertical scale factor equal
106      e3t_b (:,:,:) = e3t_n (:,:,:)
107      e3u_b (:,:,:) = e3u_n (:,:,:)
108      e3v_b (:,:,:) = e3v_n (:,:,:)
109
110      e3uw_b (:,:,:) = e3uw_n (:,:,:)
111      e3vw_b (:,:,:) = e3vw_n (:,:,:)
112      gdept_b(:,:,:) = gdept_n(:,:,:)
113      gdepw_b(:,:,:) = gdepw_n(:,:,:)
114      hu_b   (:,:)   = hu_n   (:,:)
115      hv_b   (:,:)   = hv_n   (:,:)
116      r1_hu_b(:,:)   = r1_hu_n(:,:)
117      r1_hv_b(:,:)   = r1_hv_n(:,:)
118      !
119   END SUBROUTINE iscpl_stp
120
121
122   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
123      !!----------------------------------------------------------------------
124      !!                   ***  ROUTINE iscpl_rst_interpol  ***
125      !!
126      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
127      !!                compute 2d fields of heat, salt and volume correction
128      !!
129      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
130      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
131      !!----------------------------------------------------------------------
132      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
133      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
134      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
135      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
136      !!
137      INTEGER :: ji, jj, jk, iz          !! loop index
138      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
139      !!
140      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
141      REAL(wp):: zdz, zdzm1, zdzp1
142      !!
143      REAL(wp), DIMENSION(jpi,jpj) :: zdmask , zdsmask, zvcorr, zucorr, zde3t
144      REAL(wp), DIMENSION(jpi,jpj) :: zbub   , zbvb   , zbun  , zbvn
145      REAL(wp), DIMENSION(jpi,jpj) :: zssh0  , zssh1, zhu1, zhv1
146      REAL(wp), DIMENSION(jpi,jpj) :: zsmask0, zsmask1
147      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, ztmask1, ztrp
148      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwmaskn, zwmaskb, ztmp3d
149      REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: zts0
150      !!----------------------------------------------------------------------
151
152      !! allocate variables
153
154      !! mask value to be sure
155      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
156      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
157     
158      ! compute wmask
159      zwmaskn(:,:,1) = tmask   (:,:,1)
160      zwmaskb(:,:,1) = ptmask_b(:,:,1)
161      DO jk = 2,jpk
162         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
163         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
164      END DO
165           
166      ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
167      sshb (:,:)=sshn(:,:)
168      zssh0(:,:)=sshn(:,:)
169      zsmask0(:,:) = psmask_b(:,:)
170      zsmask1(:,:) = psmask_b(:,:)
171      DO iz = 1,10    ! need to be tuned (configuration dependent) (OK for ISOMIP+)
172         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
173         DO jj = 2,jpj-1
174            DO ji = fs_2, fs_jpim1   ! vector opt.
175               jip1=ji+1; jim1=ji-1;
176               jjp1=jj+1; jjm1=jj-1;
177               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
178               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
179                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
180                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
181                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
182                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
183                  zsmask1(ji,jj)=1._wp
184               END IF
185            END DO
186         END DO
187         CALL lbc_lnk(sshn,'T',1._wp)
188         CALL lbc_lnk(zsmask1,'T',1._wp)
189         zssh0   = sshn
190         zsmask0 = zsmask1
191      END DO
192      sshn(:,:) = sshn(:,:) * ssmask(:,:)
193
194!=============================================================================
195!PM: Is this needed since introduction of VVL by default?
196      IF (.NOT.ln_linssh) THEN
197      ! Reconstruction of all vertical scale factors at now time steps
198      ! =============================================================================
199      ! Horizontal scale factor interpolations
200      ! --------------------------------------
201         DO jk = 1,jpk
202            DO jj=1,jpj
203               DO ji=1,jpi
204                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
205                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
206                  ENDIF
207               END DO
208            END DO
209         END DO
210
211         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
212         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
213         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
214
215      ! Vertical scale factor interpolations
216      ! ------------------------------------
217         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
218         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
219         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
220
221      ! t- and w- points depth
222      ! ----------------------
223         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
224         gdepw_n(:,:,1) = 0.0_wp
225         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
226         DO jj = 1,jpj
227            DO ji = 1,jpi
228               DO jk = 2,mikt(ji,jj)-1
229                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
230                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
231                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
232               END DO
233               IF (mikt(ji,jj) > 1) THEN
234                  jk = mikt(ji,jj)
235                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
236                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
237                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
238               END IF
239               DO jk = mikt(ji,jj)+1, jpk
240                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
241                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
242                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
243               END DO
244            END DO
245         END DO
246
247      ! t-, u- and v- water column thickness
248      ! ------------------------------------
249         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
250         DO jk = 1, jpkm1
251            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
252            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
253            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
254         END DO
255         !                                        ! Inverse of the local depth
256         r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
257         r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
258
259      END IF
260
261!=============================================================================
262! compute velocity
263! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
264      ub(:,:,:)=un(:,:,:)
265      vb(:,:,:)=vn(:,:,:)
266      DO jk = 1,jpk
267         DO jj = 1,jpj
268            DO ji = 1,jpi
269               un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk);
270               vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk);
271            END DO
272         END DO
273      END DO
274
275! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
276! compute barotropic velocity now and after
277      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
278      zbub(:,:)   = SUM(ztrp,DIM=3)
279      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
280      zbvb(:,:)   = SUM(ztrp,DIM=3)
281      ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 
282      zbun(:,:)   = SUM(ztrp,DIM=3)
283      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
284      zbvn(:,:)   = SUM(ztrp,DIM=3)
285
286      ! new water column
287      zhu1=0.0_wp ;
288      zhv1=0.0_wp ;
289      DO jk  = 1,jpk
290        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
291        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
292      END DO
293     
294      ! compute correction     
295      zucorr = 0._wp
296      zvcorr = 0._wp
297      DO jj = 1,jpj
298         DO ji = 1,jpi
299            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
300               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
301            END IF
302            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
303               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
304            END IF
305         END DO
306      END DO 
307     
308      ! update velocity
309      DO jk  = 1,jpk
310         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
311         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
312      END DO
313
314!=============================================================================
315      ! compute temp and salt
316      ! compute new tn and sn if we open a new cell
317      tsb (:,:,:,:) = tsn(:,:,:,:)
318      zts0(:,:,:,:) = tsn(:,:,:,:)
319      ztmask1(:,:,:) = ptmask_b(:,:,:)
320      ztmask0(:,:,:) = ptmask_b(:,:,:)
321      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
322          DO jk = 1,jpk-1
323              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
324              DO jj = 2,jpj-1
325                 DO ji = fs_2,fs_jpim1
326                      jip1=ji+1; jim1=ji-1;
327                      jjp1=jj+1; jjm1=jj-1;
328                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
329                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
330                      !! horizontal basic extrapolation
331                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
332                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
333                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
334                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
335                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
336                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
337                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
338                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
339                         ztmask1(ji,jj,jk)=1
340                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
341                      !! vertical extrapolation if horizontal extrapolation failed
342                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
343                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
344                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
345                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
346                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
347                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
348                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
349                            ztmask1(ji,jj,jk)=1._wp
350                         END IF
351                      END IF
352                  END DO
353              END DO
354          END DO
355         
356          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
357          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
358          CALL lbc_lnk(ztmask1,     'T',1._wp)
359
360          ! update
361          zts0(:,:,:,:) = tsn(:,:,:,:)
362          ztmask0 = ztmask1
363
364      END DO
365
366      ! mask new tsn field
367      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
368      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
369
370      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
371      !PM: Is this IF needed since change to VVL by default
372      IF (.NOT.ln_linssh) THEN
373         DO jk = 2,jpk-1
374            DO jj = 1,jpj
375               DO ji = 1,jpi
376                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
377                     !compute weight
378                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
379                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
380                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  ))
381                     IF (zdz .LT. 0._wp) THEN
382                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
383                     END IF
384                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
385                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
386                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk)
387                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
388                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
389                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk)
390                  END IF
391               END DO
392            END DO
393         END DO               
394      END IF
395
396      ! closed pool
397      ! -----------------------------------------------------------------------------------------
398      ! case we open a cell but no neigbour cells available to get an estimate of T and S
399      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
400         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
401         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
402         umask(:,:,:) = 0._wp
403         vmask(:,:,:) = 0._wp
404      END WHERE
405     
406      ! set mbkt and mikt to 1 in thiese location
407      WHERE (SUM(tmask,dim=3) == 0)
408         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
409         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
410      END WHERE
411      ! -------------------------------------------------------------------------------------------
412      ! compute new tn and sn if we close cell
413      ! nothing to do
414      !
415      ! deallocation tmp arrays
416      !
417   END SUBROUTINE iscpl_rst_interpol
418
419   !!======================================================================
420END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.