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/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 @ 6225

Last change on this file since 6225 was 6225, checked in by jamesharle, 8 years ago

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