source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90 @ 10001

Last change on this file since 10001 was 10001, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branch - step I.1 and I.2 (see wiki page)

File size: 20.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 phycst          ! physical constants
17   USE sbc_oce         ! surface boundary condition variables
18   USE iscplini        ! ice sheet coupling: initialisation
19   USE iscplhsb        ! ice sheet coupling: conservation
20   !
21   USE in_out_manager  ! I/O manager
22   USE iom             ! I/O module
23   USE lib_mpp         ! MPP library
24   USE lib_fortran     ! MPP library
25   USE lbclnk          ! communication
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/OCE 4.0 , NEMO Consortium (2018)
36   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
37   !! Software governed by the CeCILL licence     (./LICENSE)
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      !                       ! get restart variable
59      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
60      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
61      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
62      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
63      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:), ldxios = lrxios )  ! need to compute temperature correction
64      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b(:,:,:), ldxios = lrxios )  ! need to correct barotropic velocity
65      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b(:,:,:), ldxios = lrxios )  ! need to correct barotropic velocity
66      CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
67      !
68      CALL iscpl_init()       ! read namelist
69      !                       ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
70      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
71      !
72      IF ( ln_hsb ) THEN      ! compute correction if conservation needed
73         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
74         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
75      END IF
76         
77      !                       ! create  a domain file
78      IF( ln_meshmask .AND. ln_iscpl )   CALL dom_wri
79      !
80      IF ( ln_hsb ) THEN
81         cfile='correction'
82         cfile = TRIM( cfile )
83         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
84         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
85         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
86         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
87         CALL iom_close ( inum0 )
88      END IF
89      !
90      l_1st_euler = .TRUE.    ! next step is an euler time step
91      !
92      !                       ! set _b and _n variables equal
93      tsb (:,:,:,:) = tsn (:,:,:,:)
94      ub  (:,:,:)   = un  (:,:,:)
95      vb  (:,:,:)   = vn  (:,:,:)
96      sshb(:,:)     = sshn(:,:)
97      !
98      !                       ! set _b and _n vertical scale factor equal
99      e3t_b (:,:,:) = e3t_n (:,:,:)
100      e3u_b (:,:,:) = e3u_n (:,:,:)
101      e3v_b (:,:,:) = e3v_n (:,:,:)
102      !
103      e3uw_b (:,:,:) = e3uw_n (:,:,:)
104      e3vw_b (:,:,:) = e3vw_n (:,:,:)
105      gdept_b(:,:,:) = gdept_n(:,:,:)
106      gdepw_b(:,:,:) = gdepw_n(:,:,:)
107      hu_b   (:,:)   = hu_n   (:,:)
108      hv_b   (:,:)   = hv_n   (:,:)
109      r1_hu_b(:,:)   = r1_hu_n(:,:)
110      r1_hv_b(:,:)   = r1_hv_n(:,:)
111      !
112   END SUBROUTINE iscpl_stp
113
114
115   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
116      !!----------------------------------------------------------------------
117      !!                   ***  ROUTINE iscpl_rst_interpol  ***
118      !!
119      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
120      !!                compute 2d fields of heat, salt and volume correction
121      !!
122      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
123      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
124      !!----------------------------------------------------------------------
125      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
126      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
127      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
128      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
129      !!
130      INTEGER :: ji, jj, jk, iz          !! loop index
131      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
132      !!
133      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
134      REAL(wp):: zdz, zdzm1, zdzp1
135      !!
136      REAL(wp), DIMENSION(jpi,jpj)          :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t
137      REAL(wp), DIMENSION(jpi,jpj)          :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1
138      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn
139      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
140      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
141      REAL(wp), DIMENSION(jpi,jpj) ::   z_ssh_h0, zsshu, zsshv, zsshf
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               summsk = zsmask0(ji+1,jj)+zsmask0(ji-1,jj)+zsmask0(ji,jj+1)+zsmask0(ji,jj-1)
166               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
167                  sshn(ji,jj)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj)     &
168                  &           + zssh0(ji-1,jj)*zsmask0(ji-1,jj)     &
169                  &           + zssh0(ji,jj+1)*zsmask0(ji,jj+1)     &
170                  &           + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk
171                  zsmask1(ji,jj) = 1._wp
172               ENDIF
173            END DO
174         END DO
175         CALL lbc_lnk_multi( sshn, 'T', 1., zsmask1, 'T', 1. )
176         zssh0   = sshn
177         zsmask0 = zsmask1
178      END DO
179      sshn(:,:) = sshn(:,:) * ssmask(:,:)
180
181!=============================================================================
182!PM: Is this needed since introduction of VVL by default?
183      IF ( .NOT.ln_linssh ) THEN
184      ! Reconstruction of all vertical scale factors at now time steps
185      ! ======================================================================
186     
187!!gm Question : bug ????
188      ! at this stage  is ht_0, r1_ht0 already computed ????
189      ! I have the feeling that this should be done in a different manner...
190      ! that is iscpl_rst_interpol  should be call directly in restart !!!
191     
192      ! in my changes here I assume that ht_0 AND  r1_ht_0  have been updated
193      ! Note that the former calculation were using ht_0  so if it as not been updated ===>>> BUG
194!!gm
195     
196     
197      ! Horizontal scale factor interpolations
198      ! --------------------------------------
199         DO jj = 1, jpj
200            DO ji = 1, jpi
201               IF ( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp ) THEN
202                  DO jk = 1, jpk
203                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
204                  END DO
205               ENDIF
206            END DO
207         END DO
208         !
209!!gm  Note that if this routine is called in dom_vvl_init then all the lines below are uselss !!!
210!!        they are a duplication of dom_vvl_init lines
211
212         !                                   !==  now fields  ==!
213         !
214         !                                            !* ssh at u- and v-points)
215         DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 due to f-point
216            zsshu(ji,jj) = 0.5_wp  * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj)
217            zsshv(ji,jj) = 0.5_wp  * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj)
218            zsshf(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   & 
219               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj)
220         END DO             ;   END DO     
221         CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp , zsshf(:,:),'F', 1._wp )
222         !
223         !                                            !* hu and hv (and their inverse)
224         ht_n   (:,:) = ht_0(:,:) +  sshn(:,:)
225         hu_n   (:,:) = hu_0(:,:) + zsshu(:,:)
226         hv_n   (:,:) = hv_0(:,:) + zsshv(:,:)
227         r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )    ! ss mask mask due to ISF
228         r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
229         !
230         !                                            !* e3u, e3uw  and  e3v, e3vw
231         z_ssh_h0(:,:) = sshn (:,:) * r1_ht_0(:,:)           ! t-point
232         DO jk = 1, jpkm1
233            e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * tmask(:,:,jk) )
234         END DO
235         z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:)           ! u-point
236         DO jk = 1, jpkm1
237            e3u_n (:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) *  umask(:,:,jk) )
238            e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wumask(:,:,jk) )
239         END DO
240         z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:)           ! v-point
241         DO jk = 1, jpkm1
242            e3v_n (:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) *  vmask(:,:,jk) )
243            e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wvmask(:,:,jk) )
244         END DO
245         z_ssh_h0(:,:) = zsshf(:,:) * r1_hf_0(:,:)           ! f-point
246         DO jk = 1, jpkm1
247            e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * fmask(:,:,jk) )
248         END DO
249
250         z_ssh_h0(:,:) = 1._wp + sshn(:,:) * r1_ht_0(:,:)    ! t-point
251         !
252         IF( ln_isfcav ) THEN    ! iceshelf cavities : ssh scaling not applied over the iceshelf thickness
253            DO jk = 1, jpkm1
254               gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:)
255               gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:)
256               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
257            END DO
258         ELSE
259            DO jk = 1, jpkm1
260               gdept_n(:,:,jk) = gdept_0(:,:,jk) * z_ssh_h0(:,:)
261               gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * z_ssh_h0(:,:)
262               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:)
263            END DO
264         ENDIF
265         
266      ENDIF
267
268!=============================================================================
269! compute velocity
270! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
271      ub(:,:,:) = un(:,:,:)
272      vb(:,:,:) = vn(:,:,:)
273      DO jk = 1, jpk
274         DO jj = 1, jpj
275            DO ji = 1, jpi
276               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);
277               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);
278            END DO
279         END DO
280      END DO
281
282! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
283! compute barotropic velocity now and after
284      zbub(:,:)   = SUM( ub(:,:,:)*pe3u_b(:,:,:) , DIM=3 )
285      zbvb(:,:)   = SUM( vb(:,:,:)*pe3v_b(:,:,:) , DIM=3 )
286      zbun(:,:)   = SUM( un(:,:,:)* e3u_n(:,:,:) , DIM=3 )
287      zbvn(:,:)   = SUM( vn(:,:,:)* e3v_n(:,:,:) , DIM=3 )
288
289      ! new water column
290      zhu1=0.0_wp ;
291      zhv1=0.0_wp ;
292      DO jk  = 1,jpk
293        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
294        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
295      END DO
296     
297      ! compute correction     
298      zucorr = 0._wp
299      zvcorr = 0._wp
300      DO jj = 1, jpj
301         DO ji = 1, jpi
302            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
303               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
304            END IF
305            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
306               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
307            END IF
308         END DO
309      END DO 
310     
311      ! update velocity
312      DO jk  = 1, jpk
313         un(:,:,jk) = ( un(:,:,jk) - zucorr(:,:) ) * umask(:,:,jk)
314         vn(:,:,jk) = ( vn(:,:,jk) - zvcorr(:,:) ) * vmask(:,:,jk)
315      END DO
316
317!=============================================================================
318      ! compute temp and salt
319      ! compute new tn and sn if we open a new cell
320      tsb  (:,:,:,:) = tsn   (:,:,:,:)
321      zts0 (:,:,:,:) = tsn   (:,:,:,:)
322      ztmask1(:,:,:) = ptmask_b(:,:,:)
323      ztmask0(:,:,:) = ptmask_b(:,:,:)
324      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
325         DO jk = 1,jpk-1
326            zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
327            DO jj = 2,jpj-1
328               DO ji = fs_2,fs_jpim1
329                  jip1=ji+1; jim1=ji-1;
330                  jjp1=jj+1; jjm1=jj-1;
331                  summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
332                  IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
333                     !! horizontal basic extrapolation
334                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
335                         &            +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
336                         &            +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
337                         &            +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
338                     tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
339                         &            +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
340                         &            +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
341                         &            +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
342                     ztmask1(ji,jj,jk)=1
343                  ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
344                     !! vertical extrapolation if horizontal extrapolation failed
345                     jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
346                     summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
347                     IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
348                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
349                            &            +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / summsk
350                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
351                            &            +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / summsk
352                        ztmask1(ji,jj,jk)=1._wp
353                     ENDIF
354                  ENDIF
355               END DO
356            END DO
357         END DO
358         !
359         CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
360         !
361         ! update
362         zts0(:,:,:,:) = tsn(:,:,:,:)
363         ztmask0 = ztmask1
364         !
365      END DO
366
367      ! mask new tsn field
368      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
369      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
370
371      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
372      !PM: Is this IF needed since change to VVL by default
373      IF (.NOT.ln_linssh) THEN
374         DO jk = 2,jpk-1
375            DO jj = 1,jpj
376               DO ji = 1,jpi
377                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND.   &
378                  &      (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
379                     !compute weight
380                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
381                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
382                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  ))
383                     IF (zdz .LT. 0._wp) THEN
384                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
385                     END IF
386                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
387                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
388                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk)
389                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
390                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
391                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk)
392                  END IF
393               END DO
394            END DO
395         END DO               
396      END IF
397
398      ! closed pool
399      ! -----------------------------------------------------------------------------------------
400      ! case we open a cell but no neigbour cells available to get an estimate of T and S
401      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
402         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
403         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
404         umask(:,:,:) = 0._wp
405         vmask(:,:,:) = 0._wp
406      END WHERE
407     
408      ! set mbkt and mikt to 1 in thiese location
409      WHERE (SUM(tmask,dim=3) == 0)
410         mbkt(:,:) = 1   ;   mbku(:,:)=1   ;   mbkv(:,:) = 1
411         mikt(:,:) = 1   ;   miku(:,:)=1   ;   mikv(:,:) = 1
412      END WHERE
413      ! -------------------------------------------------------------------------------------------
414      ! compute new tn and sn if we close cell
415      ! nothing to do
416      !
417   END SUBROUTINE iscpl_rst_interpol
418
419   !!======================================================================
420END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.