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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch
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/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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   ) ! need to extrapolate T/S
61      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b   ) ! need to correct barotropic velocity
62      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b   ) ! need to correct barotropic velocity
63      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b   ) ! need to correct barotropic velocity
64      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:) )  ! need to compute temperature correction
65      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b(:,:,:) )  ! need to correct barotropic velocity
66      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b(:,:,:) )  ! need to correct barotropic velocity
67      CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:) ) ! 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( nn_msh /= 0 .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( sshn   , 'T', 1._wp )
178         CALL lbc_lnk( zsmask1, 'T', 1._wp )
179         zssh0   = sshn
180         zsmask0 = zsmask1
181      END DO
182      sshn(:,:) = sshn(:,:) * ssmask(:,:)
183
184!=============================================================================
185!PM: Is this needed since introduction of VVL by default?
186      IF ( .NOT.ln_linssh ) THEN
187      ! Reconstruction of all vertical scale factors at now time steps
188      ! =============================================================================
189      ! Horizontal scale factor interpolations
190      ! --------------------------------------
191         DO jk = 1,jpk
192            DO jj=1,jpj
193               DO ji=1,jpi
194                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
195                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) /   &
196                     &   ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
197                  ENDIF
198               END DO
199            END DO
200         END DO
201         !
202         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
203         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
204         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
205
206         ! Vertical scale factor interpolations
207         ! ------------------------------------
208         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
209         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
210         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
211         
212         ! t- and w- points depth
213         ! ----------------------
214         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
215         gdepw_n(:,:,1) = 0.0_wp
216         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
217         DO jj = 1,jpj
218            DO ji = 1,jpi
219               DO jk = 2,mikt(ji,jj)-1
220                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
221                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
222                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
223               END DO
224               IF (mikt(ji,jj) > 1) THEN
225                  jk = mikt(ji,jj)
226                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
227                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
228                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
229               END IF
230               DO jk = mikt(ji,jj)+1, jpk
231                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
232                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
233                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
234               END DO
235            END DO
236         END DO
237
238      ! t-, u- and v- water column thickness
239      ! ------------------------------------
240         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
241         DO jk = 1, jpkm1
242            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
243            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
244            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
245         END DO
246         !                                        ! Inverse of the local depth
247         r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
248         r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
249
250      END IF
251
252!=============================================================================
253! compute velocity
254! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
255      ub(:,:,:)=un(:,:,:)
256      vb(:,:,:)=vn(:,:,:)
257      DO jk = 1,jpk
258         DO jj = 1,jpj
259            DO ji = 1,jpi
260               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);
261               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);
262            END DO
263         END DO
264      END DO
265
266! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
267! compute barotropic velocity now and after
268      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
269      zbub(:,:)   = SUM(ztrp,DIM=3)
270      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
271      zbvb(:,:)   = SUM(ztrp,DIM=3)
272      ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 
273      zbun(:,:)   = SUM(ztrp,DIM=3)
274      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
275      zbvn(:,:)   = SUM(ztrp,DIM=3)
276
277      ! new water column
278      zhu1=0.0_wp ;
279      zhv1=0.0_wp ;
280      DO jk  = 1,jpk
281        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
282        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
283      END DO
284     
285      ! compute correction     
286      zucorr = 0._wp
287      zvcorr = 0._wp
288      DO jj = 1,jpj
289         DO ji = 1,jpi
290            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
291               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
292            END IF
293            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
294               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
295            END IF
296         END DO
297      END DO 
298     
299      ! update velocity
300      DO jk  = 1,jpk
301         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
302         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
303      END DO
304
305!=============================================================================
306      ! compute temp and salt
307      ! compute new tn and sn if we open a new cell
308      tsb (:,:,:,:) = tsn(:,:,:,:)
309      zts0(:,:,:,:) = tsn(:,:,:,:)
310      ztmask1(:,:,:) = ptmask_b(:,:,:)
311      ztmask0(:,:,:) = ptmask_b(:,:,:)
312      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
313          DO jk = 1,jpk-1
314              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
315              DO jj = 2,jpj-1
316                 DO ji = fs_2,fs_jpim1
317                      jip1=ji+1; jim1=ji-1;
318                      jjp1=jj+1; jjm1=jj-1;
319                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
320                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
321                      !! horizontal basic extrapolation
322                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
323                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
324                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
325                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
326                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
327                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
328                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
329                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
330                         ztmask1(ji,jj,jk)=1
331                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
332                      !! vertical extrapolation if horizontal extrapolation failed
333                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
334                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
335                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
336                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
337                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
338                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
339                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
340                            ztmask1(ji,jj,jk)=1._wp
341                         END IF
342                      END IF
343                  END DO
344              END DO
345          END DO
346         
347          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
348          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
349          CALL lbc_lnk(ztmask1,     'T',1._wp)
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.