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

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

ice sheet coupling : add comments, rm PRINT statements, 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_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          ! communication
24   USE iscplini        ! ice sheet coupling: initialisation
25   USE iscplhsb        ! ice sheet coupling: conservation
26
27   IMPLICIT NONE
28   PRIVATE
29   
30   PUBLIC   iscpl_stp          ! step management
31   PUBLIC   iscpl_rst_interpol ! routine to wet and dry
32   !!
33   !! * Substitutions 
34#  include "domzgr_substitute.h90" 
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE iscpl_stp
43      !!----------------------------------------------------------------------
44      !!                   ***  ROUTINE iscpl_stp  ***
45      !!
46      !! ** Purpose : compute initialisation
47      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
48      !!              compute correction term if needed
49      !!
50      !!----------------------------------------------------------------------
51      REAL(wp), DIMENSION(:,:  ), POINTER :: zsmask_b
52      REAL(wp), DIMENSION(:,:,:), POINTER :: ztmask_b, zumask_b, zvmask_b
53      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3t_b  , ze3u_b  , ze3v_b 
54      REAL(wp), DIMENSION(:,:,:), POINTER :: zdepw_b
55      !!----------------------------------------------------------------------
56      INTEGER  ::   inum0
57      CHARACTER(20) :: cfile
58
59      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before
60      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! e3   before
61      CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
62      CALL wrk_alloc(jpi,jpj,     zsmask_b                    )
63
64
65      !! get restart variable
66      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b   ) ! need to extrapolate T/S
67      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b   ) ! need to correct barotropic velocity
68      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b   ) ! need to correct barotropic velocity
69      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b   ) ! need to correct barotropic velocity
70      CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', ze3t_b(:,:,:) ) ! need to compute temperature correction
71      CALL iom_get( numror, jpdom_autoglo, 'fse3u_n', ze3u_b(:,:,:) ) ! need to compute volume      correction  ????
72      CALL iom_get( numror, jpdom_autoglo, 'fse3v_n', ze3v_b(:,:,:) ) ! need to compute volume      correction  ????
73      CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to compute volume      correction  ????
74
75      !! read namelist
76      CALL iscpl_init()
77
78      !!  ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
79      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
80
81      !! compute correction if conservation needed
82      IF ( ln_hsb ) THEN
83         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
84         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
85      END IF
86         
87      !! print mesh/mask
88      IF( nmsh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
89
90      IF ( ln_hsb ) THEN
91         cfile='correction'
92         cfile = TRIM( cfile )
93         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
94         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
95         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
96         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
97         CALL iom_close ( inum0 )
98      END IF
99
100      CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b ) 
101      CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b  ,ze3u_b  ,ze3v_b   ) 
102      CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b                    )
103      CALL wrk_dealloc(jpi,jpj,     zsmask_b                   )
104
105      !! next step is an euler time step
106      neuler = 0
107      !! set _b and _n variables equal
108      tsb (:,:,:,:) = tsn (:,:,:,:)
109      ub  (:,:,:  ) = un  (:,:,:  )
110      vb  (:,:,:  ) = vn  (:,:,:  )
111      sshb(:,:    ) = sshn(:,:)
112      !! set _b and _n vertical scale factor equal
113      fse3t_b (:,:,:) = fse3t_n (:,:,:)
114      fse3u_b (:,:,:) = fse3u_n (:,:,:)
115      fse3v_b (:,:,:) = fse3v_n (:,:,:)
116      IF ( lk_vvl ) THEN
117         fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
118         fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
119         fsdept_b(:,:,:) = fsdept_n(:,:,:)
120         fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
121         hu_b (:,:) = hu(:,:)
122         hv_b (:,:) = hv(:,:)
123         hur_b(:,:) = hur(:,:)
124         hvr_b(:,:) = hvr(:,:)
125      END IF
126      !
127   END SUBROUTINE iscpl_stp
128   
129   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
130      !!----------------------------------------------------------------------
131      !!                   ***  ROUTINE iscpl_rst_interpol  ***
132      !!
133      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
134      !!                compute 2d fields of heat, salt and volume correction
135      !!
136      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
137      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
138      !!----------------------------------------------------------------------
139      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
140      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
141      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
142      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
143      !!
144      INTEGER :: ji, jj, jk, iz          !! loop index
145      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
146      !!
147      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
148      REAL(wp):: zdz, zdzm1, zdzp1
149      !!
150      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
151      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
152      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, hu1, hv1
153      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
154      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
155      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
156      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
157
158      !! allocate variables
159      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
160      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d         ) 
161      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
162      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
163      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
164      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
165      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, hu1, hv1               ) 
166
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) (OK for ISOMIP+)
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._wp .AND. summsk .NE. 0._wp) 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._wp
197               END IF
198            END DO
199         END DO
200         CALL lbc_lnk(sshn,'T',1._wp)
201         CALL lbc_lnk(zsmask1,'T',1._wp)
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      ! t-, u- and v- water column thickness
260      ! ------------------------------------
261         ht(:,:) = 0._wp ; hu(:,:) = 0._wp ; hv(:,:) = 0._wp
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! compute velocity
275! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
276      ub(:,:,:)=un(:,:,:)
277      vb(:,:,:)=vn(:,:,:)
278      DO jk = 1,jpk
279         DO ji = 1,jpi
280            DO jj = 1,jpj
281               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);
282               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);
283            END DO
284         END DO
285      END DO
286
287! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
288! compute barotropic velocity now and after
289      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
290      zbub(:,:)   = SUM(ztrp,DIM=3)
291      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
292      zbvb(:,:)   = SUM(ztrp,DIM=3)
293      ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:); 
294      zbun(:,:)   = SUM(ztrp,DIM=3)
295      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
296      zbvn(:,:)   = SUM(ztrp,DIM=3)
297
298      ! new water column
299      hu1=0.0_wp ;
300      hv1=0.0_wp ;
301      DO jk  = 1,jpk
302        hu1(:,:) = hu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
303        hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
304      END DO
305     
306      ! compute correction     
307      zucorr = 0._wp
308      zvcorr = 0._wp
309      DO ji = 1,jpi
310         DO jj = 1,jpj
311            IF (zbun(ji,jj) .NE. zbub(ji,jj) .AND. hu1(ji,jj) .NE. 0._wp ) THEN
312               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/hu1(ji,jj)
313            END IF
314            IF (zbvn(ji,jj) .NE. zbvb(ji,jj) .AND. hv1(ji,jj) .NE. 0._wp ) THEN
315               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/hv1(ji,jj)
316            END IF
317         END DO
318      END DO 
319     
320      ! update velocity
321      DO jk  = 1,jpk
322         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
323         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
324      END DO
325
326!=============================================================================
327      ! compute temp and salt
328      ! compute new tn and sn if we open a new cell
329      tsb (:,:,:,:) = tsn(:,:,:,:)
330      zts0(:,:,:,:) = tsn(:,:,:,:)
331      ztmask1(:,:,:) = ptmask_b(:,:,:)
332      ztmask0(:,:,:) = ptmask_b(:,:,:)
333      DO iz = 1,10 ! resolution dependent (OK for ISOMIP+ case)
334          DO jk = 1,jpk-1
335              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
336              DO ji = 2,jpi-1
337                  DO jj = 2,jpj-1
338                      jip1=ji+1; jim1=ji-1;
339                      jjp1=jj+1; jjm1=jj-1;
340                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
341                      !! horizontal basic extrapolation
342                      IF (zdmask(ji,jj)==1._wp .AND. summsk .NE. 0._wp) THEN
343                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
344                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
345                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
346                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
347                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
348                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
349                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
350                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
351                         ztmask1(ji,jj,jk)=1
352                      END IF
353                      !! vertical extrapolation if horizontal extra failed
354                      IF (zdmask(ji,jj)==1._wp .AND. summsk==0._wp) THEN
355                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
356                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
357                         IF (zdmask(ji,jj)==1._wp .AND. summsk .NE. 0._wp ) THEN
358                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
359                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
360                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
361                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
362                            ztmask1(ji,jj,jk)=1._wp
363                         END IF
364                      END IF
365                  END DO
366              END DO
367          END DO
368         
369          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
370          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
371          CALL lbc_lnk(ztmask1,'T',1._wp)
372
373          ! update
374          zts0(:,:,:,:) = tsn(:,:,:,:)
375          ztmask0 = ztmask1
376
377      END DO
378
379      ! mask new tsn field
380      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
381      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
382
383      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
384      IF ( lk_vvl ) THEN
385         DO jk = 2,jpk-1
386            DO jj = 1,jpj
387               DO ji = 1,jpi
388                  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
389                     !compute weight
390                     zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
391                     zdz   =           fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
392                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - fsdepw_n(ji,jj,jk  ))
393                     IF (zdz .LT. 0._wp) THEN
394                        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)
395                        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)
396                        WRITE(numout,*) 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
397                        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)
398                        WRITE(numout,*) 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1)
399                        WRITE(numout,*) 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1)
400                        WRITE(numout,*) 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
401                        WRITE(numout,*) 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
402                        WRITE(numout,*) 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
403                        WRITE(numout,*) 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
404                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
405                     END IF
406                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
407                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
408                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/fse3t_n(ji,jj,jk)
409                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
410                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
411                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/fse3t_n(ji,jj,jk)
412                  END IF
413               END DO
414            END DO
415         END DO               
416      END IF
417
418      ! closed pool
419      ! -----------------------------------------------------------------------------------------
420      ! case we open a cell but no neigbour cells available to get an estimate of T and S
421      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
422         tsn(:,:,:,2)=  -99._wp  ! Special value for closed pool (checking purpose in output.init)
423         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
424         umask(:,:,:) = 0._wp
425         vmask(:,:,:) = 0._wp
426      END WHERE
427     
428      ! set mbkt and mikt to 1 in thiese location
429      WHERE (SUM(tmask,dim=3) == 0)
430         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
431         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
432      END WHERE
433      ! -------------------------------------------------------------------------------------------
434      ! compute new tn and sn if we close cell
435      ! nothing to do
436      !
437      ! deallocation tmp arrays
438      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
439      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
440      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d              ) 
441      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
442      CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
443      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
444      CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , hu1 , hv1            ) 
445
446   END SUBROUTINE iscpl_rst_interpol
447
448END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.