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/NEMO_4.0.4_GC_couple_fixrest/src/OCE/DOM – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_GC_couple_fixrest/src/OCE/DOM/iscplrst.F90 @ 15792

Last change on this file since 15792 was 15792, checked in by frrh, 2 years ago

Remove lbc_lnk calls for e3t* fields rendered redundant by
inclusion of lbc_lnk on rnf.
File size: 19.4 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. )
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( 'iscplrst', 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
202        ! call lbc_lnk("iscpl_rst_interpol", e3t_n, 'T', 1.0_wp)
203
204         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
205         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
206         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
207
208         ! Vertical scale factor interpolations
209         ! ------------------------------------
210         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
211         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
212         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
213         
214         ! t- and w- points depth
215         ! ----------------------
216         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
217         gdepw_n(:,:,1) = 0.0_wp
218         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
219         DO jj = 1,jpj
220            DO ji = 1,jpi
221               DO jk = 2,mikt(ji,jj)-1
222                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
223                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
224                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
225               END DO
226               IF (mikt(ji,jj) > 1) THEN
227                  jk = mikt(ji,jj)
228                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
229                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
230                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
231               END IF
232               DO jk = mikt(ji,jj)+1, jpk
233                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
234                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
235                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
236               END DO
237            END DO
238         END DO
239
240      ! t-, u- and v- water column thickness
241      ! ------------------------------------
242         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
243         DO jk = 1, jpkm1
244            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
245            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
246            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
247         END DO
248         !                                        ! Inverse of the local depth
249         r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
250         r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
251
252      END IF
253
254!=============================================================================
255! compute velocity
256! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
257      ub(:,:,:)=un(:,:,:)
258      vb(:,:,:)=vn(:,:,:)
259      DO jk = 1,jpk
260         DO jj = 1,jpj
261            DO ji = 1,jpi
262               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);
263               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);
264            END DO
265         END DO
266      END DO
267
268! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
269! compute barotropic velocity now and after
270      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
271      zbub(:,:)   = SUM(ztrp,DIM=3)
272      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
273      zbvb(:,:)   = SUM(ztrp,DIM=3)
274      ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 
275      zbun(:,:)   = SUM(ztrp,DIM=3)
276      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
277      zbvn(:,:)   = SUM(ztrp,DIM=3)
278
279      ! new water column
280      zhu1=0.0_wp ;
281      zhv1=0.0_wp ;
282      DO jk  = 1,jpk
283        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
284        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
285      END DO
286     
287      ! compute correction     
288      zucorr = 0._wp
289      zvcorr = 0._wp
290      DO jj = 1,jpj
291         DO ji = 1,jpi
292            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
293               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
294            END IF
295            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
296               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
297            END IF
298         END DO
299      END DO 
300     
301      ! update velocity
302      DO jk  = 1,jpk
303         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
304         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
305      END DO
306
307!=============================================================================
308      ! compute temp and salt
309      ! compute new tn and sn if we open a new cell
310      tsb (:,:,:,:) = tsn(:,:,:,:)
311      zts0(:,:,:,:) = tsn(:,:,:,:)
312      ztmask1(:,:,:) = ptmask_b(:,:,:)
313      ztmask0(:,:,:) = ptmask_b(:,:,:)
314      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
315          DO jk = 1,jpk-1
316              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
317              DO jj = 2,jpj-1
318                 DO ji = fs_2,fs_jpim1
319                      jip1=ji+1; jim1=ji-1;
320                      jjp1=jj+1; jjm1=jj-1;
321                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
322                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
323                      !! horizontal basic extrapolation
324                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
325                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
326                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
327                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
328                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
329                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
330                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
331                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
332                         ztmask1(ji,jj,jk)=1
333                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
334                      !! vertical extrapolation if horizontal extrapolation failed
335                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
336                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
337                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
338                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
339                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
340                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
341                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
342                            ztmask1(ji,jj,jk)=1._wp
343                         END IF
344                      END IF
345                  END DO
346              END DO
347          END DO
348         
349          CALL lbc_lnk_multi( 'iscplrst', tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
350
351          ! update
352          zts0(:,:,:,:) = tsn(:,:,:,:)
353          ztmask0 = ztmask1
354
355      END DO
356
357      ! mask new tsn field
358      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
359      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
360
361      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
362      !PM: Is this IF needed since change to VVL by default
363      IF (.NOT.ln_linssh) THEN
364         DO jk = 2,jpk-1
365            DO jj = 1,jpj
366               DO ji = 1,jpi
367                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND.   &
368                  &      (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
369                     !compute weight
370                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
371                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
372                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  ))
373                     IF (zdz .LT. 0._wp) THEN
374                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
375                     END IF
376                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
377                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
378                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk)
379                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
380                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
381                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk)
382                  END IF
383               END DO
384            END DO
385         END DO               
386      END IF
387
388      ! closed pool
389      ! -----------------------------------------------------------------------------------------
390      ! case we open a cell but no neigbour cells available to get an estimate of T and S
391      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
392         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
393         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
394         umask(:,:,:) = 0._wp
395         vmask(:,:,:) = 0._wp
396      END WHERE
397     
398      ! set mbkt and mikt to 1 in thiese location
399      WHERE (SUM(tmask,dim=3) == 0)
400         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
401         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
402      END WHERE
403      ! -------------------------------------------------------------------------------------------
404      ! compute new tn and sn if we close cell
405      ! nothing to do
406      !
407   END SUBROUTINE iscpl_rst_interpol
408
409   !!======================================================================
410END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.