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 @ 5920

Last change on this file since 5920 was 5920, checked in by mathiot, 8 years ago

ice sheet coupling: add treshold to define grounded area, remove useless va
riable, change some variable name + add some namelist parameter
File size: 22.0 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   PUBLIC   iscpl_rst_interpol ! routine to wet and dry
33   !!
34   !! * Substitutions 
35#  include "domzgr_substitute.h90" 
36#  include "vectopt_loop_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 correct barotropic velocity
74      CALL iom_get( numror, jpdom_autoglo, 'fse3v_n' , ze3v_b(:,:,:) )  ! need to correct barotropic velocity
75      CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to interpol vertical profile (vvl)
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
110      !! set _b and _n variables equal
111      tsb (:,:,:,:) = tsn (:,:,:,:)
112      ub  (:,:,:  ) = un  (:,:,:  )
113      vb  (:,:,:  ) = vn  (:,:,:  )
114      sshb(:,:    ) = sshn(:,:)
115
116      !! set _b and _n vertical scale factor equal
117      fse3t_b (:,:,:) = fse3t_n (:,:,:)
118      fse3u_b (:,:,:) = fse3u_n (:,:,:)
119      fse3v_b (:,:,:) = fse3v_n (:,:,:)
120      IF ( lk_vvl ) THEN
121         fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
122         fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
123         fsdept_b(:,:,:) = fsdept_n(:,:,:)
124         fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
125         hu_b (:,:) = hu(:,:)
126         hv_b (:,:) = hv(:,:)
127         hur_b(:,:) = hur(:,:)
128         hvr_b(:,:) = hvr(:,:)
129      END IF
130      !
131   END SUBROUTINE iscpl_stp
132   
133   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
134      !!----------------------------------------------------------------------
135      !!                   ***  ROUTINE iscpl_rst_interpol  ***
136      !!
137      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
138      !!                compute 2d fields of heat, salt and volume correction
139      !!
140      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
141      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
142      !!----------------------------------------------------------------------
143      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
144      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
145      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
146      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
147      !!
148      INTEGER :: ji, jj, jk, iz          !! loop index
149      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
150      !!
151      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
152      REAL(wp):: zdz, zdzm1, zdzp1
153      !!
154      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
155      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
156      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, hu1, hv1
157      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
158      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
159      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
160      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
161
162      !! allocate variables
163      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
164      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        ) 
165      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
166      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
167      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
168      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
169      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, hu1, hv1               ) 
170
171      !! mask value to be sure
172      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
173      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
174     
175      ! compute wmask
176      zwmaskn(:,:,1) = tmask   (:,:,1)
177      zwmaskb(:,:,1) = ptmask_b(:,:,1)
178      DO jk = 2,jpk
179         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
180         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
181      END DO
182           
183      ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
184      sshb (:,:)=sshn(:,:)
185      zssh0(:,:)=sshn(:,:)
186      zsmask0(:,:) = psmask_b(:,:)
187      zsmask1(:,:) = psmask_b(:,:)
188      DO iz = 1,10    ! need to be tuned (configuration dependent) (OK for ISOMIP+)
189         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
190         DO jj = 2,jpj-1
191            DO ji = fs_2, fs_jpim1   ! vector opt.
192               jip1=ji+1; jim1=ji-1;
193               jjp1=jj+1; jjm1=jj-1;
194               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
195               IF (zdsmask(ji,jj)==1._wp .AND. summsk .NE. 0._wp) THEN
196                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
197                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
198                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
199                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
200                  zsmask1(ji,jj)=1._wp
201               END IF
202            END DO
203         END DO
204         CALL lbc_lnk(sshn,'T',1._wp)
205         CALL lbc_lnk(zsmask1,'T',1._wp)
206         zssh0   = sshn
207         zsmask0 = zsmask1
208      END DO
209      sshn(:,:) = sshn(:,:) * ssmask(:,:)
210
211!=============================================================================
212      IF ( lk_vvl ) THEN
213      ! Reconstruction of all vertical scale factors at now time steps
214      ! =============================================================================
215      ! Horizontal scale factor interpolations
216      ! --------------------------------------
217         DO jk = 1,jpk
218            DO jj=1,jpj
219               DO ji=1,jpi
220                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
221                     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) )
222                  ENDIF
223               END DO
224            END DO
225         END DO
226
227         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
228         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
229         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
230
231      ! Vertical scale factor interpolations
232      ! ------------------------------------
233         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
234         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
235         CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
236
237      ! t- and w- points depth
238      ! ----------------------
239         fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
240         fsdepw_n(:,:,1) = 0.0_wp
241         fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
242         DO jj = 1,jpj
243            DO ji = 1,jpi
244               DO jk = 2,mikt(ji,jj)-1
245                  fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
246                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
247                  fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
248               END DO
249               IF (mikt(ji,jj) .GT. 1) THEN
250                  jk = mikt(ji,jj)
251                  fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
252                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
253                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
254               END IF
255               DO jk = mikt(ji,jj)+1, jpk
256                  fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
257                  fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
258                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
259               END DO
260            END DO
261         END DO
262
263      ! t-, u- and v- water column thickness
264      ! ------------------------------------
265         ht(:,:) = 0._wp ; hu(:,:) = 0._wp ; hv(:,:) = 0._wp
266         DO jk = 1, jpkm1
267            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
268            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
269            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
270         END DO
271         !                                        ! Inverse of the local depth
272         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
273         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
274
275      END IF
276
277!=============================================================================
278! compute velocity
279! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
280      ub(:,:,:)=un(:,:,:)
281      vb(:,:,:)=vn(:,:,:)
282      DO jk = 1,jpk
283         DO jj = 1,jpj
284            DO ji = 1,jpi
285               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);
286               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);
287            END DO
288         END DO
289      END DO
290
291! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
292! compute barotropic velocity now and after
293      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
294      zbub(:,:)   = SUM(ztrp,DIM=3)
295      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
296      zbvb(:,:)   = SUM(ztrp,DIM=3)
297      ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:); 
298      zbun(:,:)   = SUM(ztrp,DIM=3)
299      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
300      zbvn(:,:)   = SUM(ztrp,DIM=3)
301
302      ! new water column
303      hu1=0.0_wp ;
304      hv1=0.0_wp ;
305      DO jk  = 1,jpk
306        hu1(:,:) = hu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
307        hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
308      END DO
309     
310      ! compute correction     
311      zucorr = 0._wp
312      zvcorr = 0._wp
313      DO jj = 1,jpj
314         DO ji = 1,jpi
315            IF (zbun(ji,jj) .NE. zbub(ji,jj) .AND. hu1(ji,jj) .NE. 0._wp ) THEN
316               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/hu1(ji,jj)
317            END IF
318            IF (zbvn(ji,jj) .NE. zbvb(ji,jj) .AND. hv1(ji,jj) .NE. 0._wp ) THEN
319               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/hv1(ji,jj)
320            END IF
321         END DO
322      END DO 
323     
324      ! update velocity
325      DO jk  = 1,jpk
326         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
327         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
328      END DO
329
330!=============================================================================
331      ! compute temp and salt
332      ! compute new tn and sn if we open a new cell
333      tsb (:,:,:,:) = tsn(:,:,:,:)
334      zts0(:,:,:,:) = tsn(:,:,:,:)
335      ztmask1(:,:,:) = ptmask_b(:,:,:)
336      ztmask0(:,:,:) = ptmask_b(:,:,:)
337      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
338          DO jk = 1,jpk-1
339              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
340              DO jj = 2,jpj-1
341                 DO ji = fs_2,fs_jpim1
342                      jip1=ji+1; jim1=ji-1;
343                      jjp1=jj+1; jjm1=jj-1;
344                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
345                      IF (zdmask(ji,jj)==1._wp .AND. summsk .NE. 0._wp) THEN
346                      !! horizontal basic extrapolation
347                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
348                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
349                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
350                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
351                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
352                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
353                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
354                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
355                         ztmask1(ji,jj,jk)=1
356                      ELSEIF (zdmask(ji,jj)==1._wp .AND. summsk==0._wp) THEN
357                      !! vertical extrapolation if horizontal extrapolation failed
358                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
359                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
360                         IF (zdmask(ji,jj)==1._wp .AND. summsk .NE. 0._wp ) THEN
361                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
362                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
363                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
364                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
365                            ztmask1(ji,jj,jk)=1._wp
366                         END IF
367                      END IF
368                  END DO
369              END DO
370          END DO
371         
372          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
373          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
374          CALL lbc_lnk(ztmask1,'T',1._wp)
375
376          ! update
377          zts0(:,:,:,:) = tsn(:,:,:,:)
378          ztmask0 = ztmask1
379
380      END DO
381
382      ! mask new tsn field
383      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
384      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
385
386      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
387      IF ( lk_vvl ) 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,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
394                     zdz   =           fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
395                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - fsdepw_n(ji,jj,jk  ))
396                     IF (zdz .LT. 0._wp) THEN
397                        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)
398                        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)
399                        WRITE(numout,*) 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
400                        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)
401                        WRITE(numout,*) 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1)
402                        WRITE(numout,*) 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1)
403                        WRITE(numout,*) 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
404                        WRITE(numout,*) 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
405                        WRITE(numout,*) 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
406                        WRITE(numout,*) 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
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)      )/fse3t_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)      )/fse3t_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  , hu1 , hv1            ) 
448
449   END SUBROUTINE iscpl_rst_interpol
450
451END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.