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.
iscplrst.F90 in NEMO/branches/UKMO/dev_r9950_GO8_package/src/OCE/DOM – NEMO

source: NEMO/branches/UKMO/dev_r9950_GO8_package/src/OCE/DOM/iscplrst.F90 @ 10326

Last change on this file since 10326 was 10326, checked in by davestorkey, 5 years ago

UKMO/dev_r9950_GO8_package: clear SVN keywords.

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