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

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

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