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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplrst.F90 @ 10978

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

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : DOM and sshwzv.F90. (sshwzv.F90 will be rewritten somewhat when we rewrite the time-level swapping but I've done it anyway). Passes all non-AGRIF SETTE tests.
  • Property svn:keywords set to Id
File size: 19.9 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( Kbb, Kmm )
43      !!----------------------------------------------------------------------
44      !!                   ***  ROUTINE iscpl_stp  ***
45      !!
46      !! ** Purpose : compute initialisation
47      !!              compute extrapolation of restart variable uu(Kmm), vv(Kmm), ts(Kmm), ssh(Kmm) (wetting/drying)   
48      !!              compute correction term if needed
49      !!
50      !!----------------------------------------------------------------------
51      INTEGER, INTENT(in) :: Kbb, Kmm   ! time level indices
52      !
53      INTEGER  ::   inum0
54      REAL(wp), DIMENSION(jpi,jpj)     ::   zsmask_b
55      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztmask_b, zumask_b, zvmask_b
56      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t_b  , ze3u_b  , ze3v_b 
57      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepw_b
58      CHARACTER(20) :: cfile
59      !!----------------------------------------------------------------------
60      !
61      !                       ! get restart variable
62      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
63      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
64      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
65      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
66      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:), ldxios = lrxios )  ! need to compute temperature correction
67      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b(:,:,:), ldxios = lrxios )  ! need to correct barotropic velocity
68      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b(:,:,:), ldxios = lrxios )  ! need to correct barotropic velocity
69      CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
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( Kbb, Kmm, 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( Kbb, Kmm, 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      ts (:,:,:,:,Kbb) = ts (:,:,:,:,Kmm)
97      uu  (:,:,:,Kbb)   = uu  (:,:,:,Kmm)
98      vv  (:,:,:,Kbb)   = vv  (:,:,:,Kmm)
99      ssh(:,:,Kbb)     = ssh(:,:,Kmm)
100      !
101      !                       ! set _b and _n vertical scale factor equal
102      e3t (:,:,:,Kbb) = e3t (:,:,:,Kmm)
103      e3u (:,:,:,Kbb) = e3u (:,:,:,Kmm)
104      e3v (:,:,:,Kbb) = e3v (:,:,:,Kmm)
105      !
106      e3uw (:,:,:,Kbb) = e3uw (:,:,:,Kmm)
107      e3vw (:,:,:,Kbb) = e3vw (:,:,:,Kmm)
108      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
109      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
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 ( Kbb, Kmm, 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 ts(Kmm), uu(Kmm), vv(Kmm) and ssh(Kmm) in case of evolving geometry of ice shelves
123      !!                compute 2d fields of heat, salt and volume correction
124      !!
125      !! ** Method  :   ts(Kmm) : extrapolation from neigbourg cells
126      !!                uu(Kmm), vv(Kmm) : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
127      !!----------------------------------------------------------------------
128      INTEGER                     , INTENT(in ) :: Kbb, Kmm                        !! time level indices
129      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
130      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
131      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
132      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
133      !!
134      INTEGER :: ji, jj, jk, iz          !! loop index
135      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
136      !!
137      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
138      REAL(wp):: zdz, zdzm1, zdzp1
139      !!
140      REAL(wp), DIMENSION(jpi,jpj)          :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t
141      REAL(wp), DIMENSION(jpi,jpj)          :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1
142      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn, ztrp
143      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
144      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
145      !!----------------------------------------------------------------------
146      !
147      !                 ! mask value to be sure
148      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) * ptmask_b(:,:,:)
149      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) * ptmask_b(:,:,:)
150      !
151      !                 ! compute wmask
152      zwmaskn(:,:,1) = tmask   (:,:,1)
153      zwmaskb(:,:,1) = ptmask_b(:,:,1)
154      DO jk = 2,jpk
155         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
156         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
157      END DO
158      !   
159      !                 ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
160      ssh (:,:,Kbb)=ssh(:,:,Kmm)
161      zssh0(:,:)=ssh(:,:,Kmm)
162      zsmask0(:,:) = psmask_b(:,:)
163      zsmask1(:,:) = psmask_b(:,:)
164      DO iz = 1, 10                 ! need to be tuned (configuration dependent) (OK for ISOMIP+)
165         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
166         DO jj = 2,jpj-1
167            DO ji = fs_2, fs_jpim1   ! vector opt.
168               jip1=ji+1; jim1=ji-1;
169               jjp1=jj+1; jjm1=jj-1;
170               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
171               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
172                  ssh(ji,jj,Kmm)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
173                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
174                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
175                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
176                  zsmask1(ji,jj)=1._wp
177               ENDIF
178            END DO
179         END DO
180         CALL lbc_lnk_multi( 'iscplrst', ssh(:,:,Kmm), 'T', 1., zsmask1, 'T', 1. )
181         zssh0   = ssh(:,:,Kmm)
182         zsmask0 = zsmask1
183      END DO
184      ssh(:,:,Kmm) = ssh(:,:,Kmm) * ssmask(:,:)
185
186!=============================================================================
187!PM: Is this needed since introduction of VVL by default?
188      IF ( .NOT.ln_linssh ) THEN
189      ! Reconstruction of all vertical scale factors at now time steps
190      ! =============================================================================
191      ! Horizontal scale factor interpolations
192      ! --------------------------------------
193         DO jk = 1,jpk
194            DO jj=1,jpj
195               DO ji=1,jpi
196                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
197                     e3t(ji,jj,jk,Kmm) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kmm) /   &
198                     &   ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
199                  ENDIF
200               END DO
201            END DO
202         END DO
203         !
204         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
205         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
206         CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
207
208         ! Vertical scale factor interpolations
209         ! ------------------------------------
210         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
211         CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
212         CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
213         
214         ! t- and w- points depth
215         ! ----------------------
216         gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
217         gdepw(:,:,1,Kmm) = 0.0_wp
218         gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
219         DO jj = 1,jpj
220            DO ji = 1,jpi
221               DO jk = 2,mikt(ji,jj)-1
222                  gdept(ji,jj,jk,Kmm) = gdept_0(ji,jj,jk)
223                  gdepw(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk)
224                  gde3w(ji,jj,jk) = gdept_0(ji,jj,jk) - ssh(ji,jj,Kmm)
225               END DO
226               IF (mikt(ji,jj) > 1) THEN
227                  jk = mikt(ji,jj)
228                  gdept(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w(ji,jj,jk,Kmm)
229                  gdepw(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk)
230                  gde3w(ji,jj,jk) = gdept(ji,jj,jk  ,Kmm) - ssh   (ji,jj,Kmm)
231               END IF
232               DO jk = mikt(ji,jj)+1, jpk
233                  gdept(ji,jj,jk,Kmm) = gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)
234                  gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
235                  gde3w(ji,jj,jk) = gdept(ji,jj,jk  ,Kmm) - ssh (ji,jj,Kmm)
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(:,:,jk,Kmm) * umask(:,:,jk)
245            hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
246            ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * 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      uu(:,:,:,Kbb)=uu(:,:,:,Kmm)
258      vv(:,:,:,Kbb)=vv(:,:,:,Kmm)
259      DO jk = 1,jpk
260         DO jj = 1,jpj
261            DO ji = 1,jpi
262               uu(ji,jj,jk,Kmm) = uu(ji,jj,jk,Kbb)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u(ji,jj,jk,Kmm)*umask(ji,jj,jk);
263               vv(ji,jj,jk,Kmm) = vv(ji,jj,jk,Kbb)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v(ji,jj,jk,Kmm)*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(:,:,:) = uu(:,:,:,Kbb)*pe3u_b(:,:,:); 
271      zbub(:,:)   = SUM(ztrp,DIM=3)
272      ztrp(:,:,:) = vv(:,:,:,Kbb)*pe3v_b(:,:,:); 
273      zbvb(:,:)   = SUM(ztrp,DIM=3)
274      ztrp(:,:,:) = uu(:,:,:,Kmm)*e3u(:,:,:,Kmm); 
275      zbun(:,:)   = SUM(ztrp,DIM=3)
276      ztrp(:,:,:) = vv(:,:,:,Kmm)*e3v(:,:,:,Kmm); 
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(:,:,jk,Kmm) * umask(:,:,jk)
284        zhv1(:,:) = zhv1(:,:) + e3v(:,:,jk,Kmm) * 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         uu(:,:,jk,Kmm)=(uu(:,:,jk,Kmm) - zucorr(:,:))*umask(:,:,jk)
304         vv(:,:,jk,Kmm)=(vv(:,:,jk,Kmm) - 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      ts (:,:,:,:,Kbb) = ts(:,:,:,:,Kmm)
311      zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
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                         ts(ji,jj,jk,1,Kmm)=( 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                         ts(ji,jj,jk,2,Kmm)=( 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                            ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
339                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
340                            ts(ji,jj,jk,2,Kmm)=( 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', ts(:,:,:,jp_tem,Kmm), 'T', 1., ts(:,:,:,jp_sal,Kmm), 'T', 1., ztmask1, 'T', 1.)
350
351          ! update
352          zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
353          ztmask0 = ztmask1
354
355      END DO
356
357      ! mask new ts(:,:,:,:,Kmm) field
358      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) * tmask(:,:,:)
359      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) * 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(ji,jj,jk+1,Kmm) - pdepw_b(ji,jj,jk+1))
371                     zdz   =           gdepw(ji,jj,jk+1,Kmm) - pdepw_b(ji,jj,jk  ) 
372                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw(ji,jj,jk  ,Kmm))
373                     IF (zdz .LT. 0._wp) THEN
374                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
375                     END IF
376                     ts(ji,jj,jk,jp_tem,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_tem,Kbb)      &
377                        &                   + zdz  *ts(ji,jj,jk  ,jp_tem,Kbb)      &
378                        &                   + zdzm1*ts(ji,jj,jk-1,jp_tem,Kbb)      )/e3t(ji,jj,jk,Kmm)
379                     ts(ji,jj,jk,jp_sal,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_sal,Kbb)      &
380                        &                   + zdz  *ts(ji,jj,jk  ,jp_sal,Kbb)      &
381                        &                   + zdzm1*ts(ji,jj,jk-1,jp_sal,Kbb)      )/e3t(ji,jj,jk,Kmm)
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. ts(:,:,:,2,Kmm) == 0._wp) 
392         ts(:,:,:,2,Kmm) = -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.