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/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/TOOLS/DOMAINcfg/src – NEMO

source: branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/TOOLS/DOMAINcfg/src/iscplrst.f90 @ 8749

Last change on this file since 8749 was 8749, checked in by jcastill, 6 years ago

Remove svn keywords

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