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

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

ice sheet coupling branche: cosmetic changes
File size: 21.9 KB
Line 
1MODULE iscplrst
2   !!======================================================================
3   !!                       ***  MODULE  iscplrst***
4   !! Ocean forcing: update the restart file in case of ice sheet/ocean coupling
5   !!=====================================================================
6   !! History :  NEMO  ! 2015-01 P. Mathiot: original
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   iscpl_stp          : step management
11   !!   iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
12   !!----------------------------------------------------------------------
13   USE 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   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE iscpl_stp
44      !!----------------------------------------------------------------------
45      !!                   ***  ROUTINE iscpl_stp  ***
46      !!
47      !! ** Purpose : compute initialisation
48      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
49      !!              compute correction term if needed
50      !!
51      !!----------------------------------------------------------------------
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      !!----------------------------------------------------------------------
57      INTEGER  ::   inum0
58      CHARACTER(20) :: cfile
59
60      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before
61      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! e3   before
62      CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
63      CALL wrk_alloc(jpi,jpj,     zsmask_b                    )
64
65
66      !! get restart variable
67      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b   ) ! need to extrapolate T/S
68      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b   ) ! need to correct barotropic velocity
69      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b   ) ! need to correct barotropic velocity
70      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b   ) ! need to correct barotropic velocity
71      CALL iom_get( numror, jpdom_autoglo, 'fse3t_n' , ze3t_b(:,:,:) )  ! need to compute temperature correction
72      CALL iom_get( numror, jpdom_autoglo, 'fse3u_n' , ze3u_b(:,:,:) )  ! need to correct barotropic velocity
73      CALL iom_get( numror, jpdom_autoglo, 'fse3v_n' , ze3v_b(:,:,:) )  ! need to correct barotropic velocity
74      CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to interpol vertical profile (vvl)
75
76      !! read namelist
77      CALL iscpl_init()
78
79      !!  ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
80      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
81
82      !! compute correction if conservation needed
83      IF ( ln_hsb ) THEN
84         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
85         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
86      END IF
87         
88      !! print mesh/mask
89      IF( nmsh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
90
91      IF ( ln_hsb ) THEN
92         cfile='correction'
93         cfile = TRIM( cfile )
94         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
95         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
96         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
97         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
98         CALL iom_close ( inum0 )
99      END IF
100
101      CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b ) 
102      CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b  ,ze3u_b  ,ze3v_b   ) 
103      CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b                    )
104      CALL wrk_dealloc(jpi,jpj,     zsmask_b                   )
105
106      !! next step is an euler time step
107      neuler = 0
108
109      !! set _b and _n variables equal
110      tsb (:,:,:,:) = tsn (:,:,:,:)
111      ub  (:,:,:  ) = un  (:,:,:  )
112      vb  (:,:,:  ) = vn  (:,:,:  )
113      sshb(:,:    ) = sshn(:,:)
114
115      !! set _b and _n vertical scale factor equal
116      fse3t_b (:,:,:) = fse3t_n (:,:,:)
117      fse3u_b (:,:,:) = fse3u_n (:,:,:)
118      fse3v_b (:,:,:) = fse3v_n (:,:,:)
119      IF ( lk_vvl ) THEN
120         fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
121         fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
122         fsdept_b(:,:,:) = fsdept_n(:,:,:)
123         fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
124         hu_b (:,:) = hu(:,:)
125         hv_b (:,:) = hv(:,:)
126         hur_b(:,:) = hur(:,:)
127         hvr_b(:,:) = hvr(:,:)
128      END IF
129      !
130   END SUBROUTINE iscpl_stp
131   
132   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
133      !!----------------------------------------------------------------------
134      !!                   ***  ROUTINE iscpl_rst_interpol  ***
135      !!
136      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
137      !!                compute 2d fields of heat, salt and volume correction
138      !!
139      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
140      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
141      !!----------------------------------------------------------------------
142      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
143      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
144      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
145      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
146      !!
147      INTEGER :: ji, jj, jk, iz          !! loop index
148      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
149      !!
150      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
151      REAL(wp):: zdz, zdzm1, zdzp1
152      !!
153      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
154      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
155      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, hu1, hv1
156      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
157      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
158      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
159      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
160
161      !! allocate variables
162      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
163      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        ) 
164      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
165      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
166      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
167      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
168      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, hu1, hv1               ) 
169
170      !! mask value to be sure
171      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
172      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
173     
174      ! compute wmask
175      zwmaskn(:,:,1) = tmask   (:,:,1)
176      zwmaskb(:,:,1) = ptmask_b(:,:,1)
177      DO jk = 2,jpk
178         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
179         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
180      END DO
181           
182      ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
183      sshb (:,:)=sshn(:,:)
184      zssh0(:,:)=sshn(:,:)
185      zsmask0(:,:) = psmask_b(:,:)
186      zsmask1(:,:) = psmask_b(:,:)
187      DO iz = 1,10    ! need to be tuned (configuration dependent) (OK for ISOMIP+)
188         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
189         DO ji = 2,jpi-1
190            DO jj = 2,jpj-1
191               jip1=ji+1; jim1=ji-1;
192               jjp1=jj+1; jjm1=jj-1;
193               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
194               IF (zdsmask(ji,jj)==1._wp .AND. summsk .NE. 0._wp) THEN
195                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
196                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
197                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
198                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
199                  zsmask1(ji,jj)=1._wp
200               END IF
201            END DO
202         END DO
203         CALL lbc_lnk(sshn,'T',1._wp)
204         CALL lbc_lnk(zsmask1,'T',1._wp)
205         zssh0   = sshn
206         zsmask0 = zsmask1
207      END DO
208      sshn(:,:) = sshn(:,:) * ssmask(:,:)
209
210!=============================================================================
211      IF ( lk_vvl ) 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                     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) )
221                  ENDIF
222               END DO
223            END DO
224         END DO
225
226         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
227         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
228         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
229
230      ! Vertical scale factor interpolations
231      ! ------------------------------------
232         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
233         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
234         CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
235
236      ! t- and w- points depth
237      ! ----------------------
238         fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
239         fsdepw_n(:,:,1) = 0.0_wp
240         fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
241         DO jj = 1,jpj
242            DO ji = 1,jpi
243               DO jk = 2,mikt(ji,jj)-1
244                  fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
245                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
246                  fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
247               END DO
248               IF (mikt(ji,jj) .GT. 1) THEN
249                  jk = mikt(ji,jj)
250                  fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
251                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
252                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
253               END IF
254               DO jk = mikt(ji,jj)+1, jpk
255                  fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
256                  fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
257                  fsde3w_n(ji,jj,jk) = fsdept_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(:,:) = 0._wp ; hu(:,:) = 0._wp ; hv(:,:) = 0._wp
265         DO jk = 1, jpkm1
266            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
267            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
268            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
269         END DO
270         !                                        ! Inverse of the local depth
271         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
272         hvr(:,:) = 1._wp / ( hv(:,:) + 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 ji = 1,jpi
283            DO jj = 1,jpj
284               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);
285               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);
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(:,:,:)*fse3u_n(:,:,:); 
297      zbun(:,:)   = SUM(ztrp,DIM=3)
298      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
299      zbvn(:,:)   = SUM(ztrp,DIM=3)
300
301      ! new water column
302      hu1=0.0_wp ;
303      hv1=0.0_wp ;
304      DO jk  = 1,jpk
305        hu1(:,:) = hu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
306        hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
307      END DO
308     
309      ! compute correction     
310      zucorr = 0._wp
311      zvcorr = 0._wp
312      DO ji = 1,jpi
313         DO jj = 1,jpj
314            IF (zbun(ji,jj) .NE. zbub(ji,jj) .AND. hu1(ji,jj) .NE. 0._wp ) THEN
315               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/hu1(ji,jj)
316            END IF
317            IF (zbvn(ji,jj) .NE. zbvb(ji,jj) .AND. hv1(ji,jj) .NE. 0._wp ) THEN
318               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/hv1(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,10 ! resolution dependent (OK for ISOMIP+ case)
337          DO jk = 1,jpk-1
338              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
339              DO ji = 2,jpi-1
340                  DO jj = 2,jpj-1
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                      !! horizontal basic extrapolation
345                      IF (zdmask(ji,jj)==1._wp .AND. summsk .NE. 0._wp) THEN
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                      END IF
356                      !! vertical extrapolation if horizontal extra failed
357                      IF (zdmask(ji,jj)==1._wp .AND. summsk==0._wp) THEN
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.