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/2019/dev_r11351_fldread_with_XIOS/src/OCE/DOM – NEMO

source: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/DOM/iscplrst.F90 @ 11405

Last change on this file since 11405 was 11405, checked in by andmirek, 5 years ago

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