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/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 @ 5820

Last change on this file since 5820 was 5820, checked in by mathiot, 9 years ago

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