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

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

ice sheet coupling: reproducibility fixed in MISOMIP configuration + contrain bathy to be constant during all the run
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_rst_interpol : restart interpolation in case of coupling with ice sheet
11   !!----------------------------------------------------------------------
12   USE dom_oce         ! ocean space and time domain
13   USE domwri          ! ocean space and time domain
14   USE domvvl, ONLY : dom_vvl_interpol
15   USE phycst          ! physical constants
16   USE sbc_oce         ! surface boundary condition variables
17   USE oce             ! global tra/dyn variable
18   USE in_out_manager  ! I/O manager
19   USE iom             ! I/O module
20   USE lib_mpp         ! MPP library
21   USE lib_fortran     ! MPP library
22   USE wrk_nemo        ! Memory allocation
23   USE lbclnk
24   USE domngb
25   USE sbc_ice, ONLY : lk_lim3
26   USE iscplini
27   USE iscplhsb
28
29   IMPLICIT NONE
30   PRIVATE
31   
32   PUBLIC   iscpl_stp       
33   PUBLIC   iscpl_rst_interpol ! routine to wet and dry
34   !!
35   !! * Substitutions 
36#  include "domzgr_substitute.h90" 
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE iscpl_stp
45      !!----------------------------------------------------------------------
46      !!                   ***  ROUTINE iscpl_stp  ***
47      !!
48      !! ** Purpose : compute initialisation
49      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
50      !!              compute correction term if needed
51      !!
52      !!----------------------------------------------------------------------
53      REAL(wp), DIMENSION(:,:  ), POINTER :: zsmask_b
54      REAL(wp), DIMENSION(:,:,:), POINTER :: ztmask_b, zumask_b, zvmask_b
55      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3t_b  , ze3u_b  , ze3v_b 
56      REAL(wp), DIMENSION(:,:,:), POINTER :: zdepw_b
57!!----------------------------------------------------------------------
58      INTEGER  ::   inum0
59      CHARACTER(20) :: cfile
60
61      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before
62      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! e3   before
63      CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
64      CALL wrk_alloc(jpi,jpj,     zsmask_b                    )
65
66
67      !! get restart variable
68      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b   ) ! need to extrapolate T/S
69      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b   ) ! need to correct barotropic velocity
70      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b   ) ! need to correct barotropic velocity
71      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b   ) ! need to correct barotropic velocity
72      CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', ze3t_b(:,:,:) ) ! need to compute temperature correction
73      CALL iom_get( numror, jpdom_autoglo, 'fse3u_n', ze3u_b(:,:,:) ) ! need to compute volume      correction  ????
74      CALL iom_get( numror, jpdom_autoglo, 'fse3v_n', ze3v_b(:,:,:) ) ! need to compute volume      correction  ????
75      CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to compute volume      correction  ????
76
77      !! read namelist
78      CALL iscpl_init()
79
80      !!  ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
81      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
82
83      !! compute correction if conservation needed
84      IF ( ln_hsb ) THEN
85         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
86         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
87      END IF
88         
89      !! print mesh/mask
90      IF( nmsh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
91
92      IF ( ln_hsb ) THEN
93         cfile='correction'
94         cfile = TRIM( cfile )
95         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
96         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
97         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
98         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
99         CALL iom_rstput( 0, 0, inum0, 'e3t_diff', (e3t_0(:,:,:)*tmask(:,:,:) - ze3t_b(:,:,:)*ztmask_b(:,:,:))*tmask(:,:,:) )
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      !! set _b and _n variables equal
111      tsb (:,:,:,:) = tsn (:,:,:,:)
112      ub  (:,:,:  ) = un  (:,:,:  )
113      vb  (:,:,:  ) = vn  (:,:,:  )
114      sshb(:,:    ) = sshn(:,:)
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
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      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
162      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d         ) 
163      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
164      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
165      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
166      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
167      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, hu1, hv1               ) 
168      !! mask value to be sure
169      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
170      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
171     
172      ! compute wmask
173      zwmaskn(:,:,1) = tmask   (:,:,1)
174      zwmaskb(:,:,1) = ptmask_b(:,:,1)
175      DO jk = 2,jpk
176         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
177         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
178      END DO
179           
180      ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
181      sshb (:,:)=sshn(:,:)
182      zssh0(:,:)=sshn(:,:)
183      zsmask0(:,:) = psmask_b(:,:)
184      zsmask1(:,:) = psmask_b(:,:)
185      DO iz = 1,10    ! need to be tuned (configuration dependent)
186         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
187         DO ji = 2,jpi-1
188            DO jj = 2,jpj-1
189               jip1=ji+1; jim1=ji-1;
190               jjp1=jj+1; jjm1=jj-1;
191               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
192               IF (zdsmask(ji,jj)==1 .AND. summsk .NE. 0) THEN
193                  PRINT *, 'add ssh to 1 cell',ji,jj,narea
194                  sshn(ji,jj)=(zssh0(jip1,jj)*zsmask0(jip1,jj)     &
195                  &           +zssh0(jim1,jj)*zsmask0(jim1,jj)     &
196                  &           +zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
197                  &           +zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
198                  zsmask1(ji,jj)=1
199               END IF
200            END DO
201         END DO
202         CALL lbc_lnk(sshn,'T',1.)
203         CALL lbc_lnk(zsmask1,'T',1.)
204         zssh0   = sshn
205         zsmask0 = zsmask1
206      END DO
207      sshn(:,:) = sshn(:,:) * ssmask(:,:)
208
209!=============================================================================
210      IF ( lk_vvl ) THEN
211! compute fse3t_n
212         DO jk = 1,jpk
213            fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) )
214         END DO
215
216! compute fse3u/v ... (call interpolation vvl)
217      ! Reconstruction of all vertical scale factors at now and before time steps
218      ! =============================================================================
219      ! Horizontal scale factor interpolations
220      ! --------------------------------------
221         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
222         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
223         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
224
225      ! Vertical scale factor interpolations
226      ! ------------------------------------
227         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
228         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
229         CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
230
231      ! t- and w- points depth
232      ! ----------------------
233         fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
234         fsdepw_n(:,:,1) = 0.0_wp
235         fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
236         DO jj = 1,jpj
237            DO ji = 1,jpi
238               DO jk = 2,mikt(ji,jj)-1
239                  fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
240                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
241                  fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
242               END DO
243               IF (mikt(ji,jj) .GT. 1) THEN
244                  jk = mikt(ji,jj)
245                  fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
246                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
247                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
248               END IF
249               DO jk = mikt(ji,jj)+1, jpk
250                  fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
251                  fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
252                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
253               END DO
254            END DO
255         END DO
256
257         hu(:,:) = 0._wp                          ! Ocean depth at U-points
258         hv(:,:) = 0._wp                          ! Ocean depth at V-points
259         ht(:,:) = 0._wp                          ! Ocean depth at T-points
260         DO jk = 1, jpkm1
261            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
262            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
263            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
264         END DO
265         !                                        ! Inverse of the local depth
266         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
267         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
268
269      END IF
270
271!=============================================================================
272
273! compute velocity
274! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
275      ub(:,:,:)=un(:,:,:)
276      vb(:,:,:)=vn(:,:,:)
277      DO jk = 1,jpk
278         DO ji = 1,jpi
279            DO jj = 1,jpj
280               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);
281               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);
282            END DO
283         END DO
284      END DO
285
286! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
287! compute barotropic velocity now and after
288      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
289      zbub(:,:)   = SUM(ztrp,DIM=3)
290      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
291      zbvb(:,:)   = SUM(ztrp,DIM=3)
292      ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:); 
293      zbun(:,:)   = SUM(ztrp,DIM=3)
294      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
295      zbvn(:,:)   = SUM(ztrp,DIM=3)
296
297! Already know ????????
298      hu1=0.0_wp ;
299      hv1=0.0_wp ;
300      DO jk  = 1,jpk
301        hu1(:,:) = hu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
302        hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
303      END DO
304      zucorr = 0._wp
305      zvcorr = 0._wp
306      DO ji = 1,jpi
307         DO jj = 1,jpj
308            IF (zbun(ji,jj) .NE. zbub(ji,jj) .AND. hu1(ji,jj) .NE. 0._wp ) THEN
309               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/hu1(ji,jj)
310            END IF
311            IF (zbvn(ji,jj) .NE. zbvb(ji,jj) .AND. hv1(ji,jj) .NE. 0._wp ) THEN
312               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/hv1(ji,jj)
313            END IF
314         END DO
315      END DO
316      DO jk  = 1,jpk
317         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
318         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
319      END DO
320
321!=============================================================================
322      ! compute temp and salt
323      ! compute new tn and sn if we open a new cell
324      tsb (:,:,:,:) = tsn(:,:,:,:)
325      zts0(:,:,:,:) = tsn(:,:,:,:)
326      ztmask1(:,:,:) = ptmask_b(:,:,:)
327      ztmask0(:,:,:) = ptmask_b(:,:,:)
328      DO iz = 1,10
329          DO jk = 1,jpk-1
330              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
331              DO ji = 2,jpi-1
332                  DO jj = 2,jpj-1
333                      jip1=ji+1; jim1=ji-1;
334                      jjp1=jj+1; jjm1=jj-1;
335                      summsk=(ztmask0(jip1,jj,jk)+ztmask0(jim1,jj,jk)+ztmask0(ji,jjp1,jk)+ztmask0(ji,jjm1,jk))
336                      !! horizontal extrapolation
337                      IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0) THEN
338                         tsn(ji,jj,jk,1)=( zts0(jip1,jj,jk,1)*ztmask0(jip1,jj,jk)     &
339                         &                +zts0(jim1,jj,jk,1)*ztmask0(jim1,jj,jk)     &
340                         &                +zts0(ji,jjp1,jk,1)*ztmask0(ji,jjp1,jk)     &
341                         &                +zts0(ji,jjm1,jk,1)*ztmask0(ji,jjm1,jk))/summsk
342                         tsn(ji,jj,jk,2)=( zts0(jip1,jj,jk,2)*ztmask0(jip1,jj,jk)     &
343                         &                +zts0(jim1,jj,jk,2)*ztmask0(jim1,jj,jk)     &
344                         &                +zts0(ji,jjp1,jk,2)*ztmask0(ji,jjp1,jk)     &
345                         &                +zts0(ji,jjm1,jk,2)*ztmask0(ji,jjm1,jk))/summsk
346                         ztmask1(ji,jj,jk)=1
347                      END IF
348                      !! vertical extrapolation if horizontal extra failed
349                      IF (zdmask(ji,jj)==1 .AND. summsk==0) THEN
350                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
351                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
352                         IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0 ) THEN
353                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
354                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
355                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
356                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
357                            ztmask1(ji,jj,jk)=1
358                         END IF
359                      END IF
360                      !! horizontal corner extrapolation if horizontal and vertical failed
361                      IF (zdmask(ji,jj)==1 .AND. summsk==0) THEN
362                         jip1=ji+1; jim1=ji-1;
363                         jjp1=jj+1; jjm1=jj-1;
364                         summsk=(ztmask0(jip1,jjp1,jk)+ztmask0(jim1,jjm1,jk)+ztmask0(jim1,jjp1,jk)+ztmask0(jip1,jjm1,jk))
365                         IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0) THEN
366                            tsn(ji,jj,jk,1)=( zts0(jip1,jjp1,jk,1)*ztmask0(jip1,jjp1,jk)     &
367                            &                +zts0(jim1,jjm1,jk,1)*ztmask0(jim1,jjm1,jk)     &
368                            &                +zts0(jim1,jjp1,jk,1)*ztmask0(jim1,jjp1,jk)     &
369                            &                +zts0(jip1,jjm1,jk,1)*ztmask0(jip1,jjm1,jk))/summsk
370                            tsn(ji,jj,jk,2)=( zts0(jip1,jjp1,jk,2)*ztmask0(jip1,jjp1,jk)     &
371                            &                +zts0(jim1,jjm1,jk,2)*ztmask0(jim1,jjm1,jk)     &
372                            &                +zts0(jim1,jjp1,jk,2)*ztmask0(jim1,jjp1,jk)     &
373                            &                +zts0(jip1,jjm1,jk,2)*ztmask0(jip1,jjm1,jk))/summsk
374                            ztmask1(ji,jj,jk)=1
375                         END IF
376                      END IF
377                  END DO
378              END DO
379          END DO
380          CALL lbc_lnk(tsn(:,:,:,1),'T',1.)
381          CALL lbc_lnk(tsn(:,:,:,2),'T',1.)
382          CALL lbc_lnk(ztmask1,'T',1.)
383          zts0(:,:,:,:) = tsn(:,:,:,:)
384          ztmask0 = ztmask1
385      END DO
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.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.0_wp) THEN
400                        PRINT *, 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1)
401                        PRINT *, 'ERROR dz n             = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj)
402                        PRINT *, 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
403                        PRINT *, 'ERROR dz b             = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj)
404                        PRINT *, 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1)
405                        PRINT *, 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1)
406                        PRINT *, 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
407                        PRINT *, 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
408                        PRINT *, 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
409                        PRINT *, 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
410                        STOP
411                     END IF
412                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
413                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
414                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/fse3t_n(ji,jj,jk)
415                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
416                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
417                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/fse3t_n(ji,jj,jk)
418                  END IF
419               END DO
420            END DO
421         END DO               
422      END IF
423
424      ! Special value for closed pool and set the mask to 0 to run
425      WHERE (tmask(:,:,:) == 1.0 .AND. tsn(:,:,:,2) == 0._wp) 
426         tsn(:,:,:,2)=  -99._wp
427         tmask(:,:,:) = 0.0
428         umask(:,:,:) = 0.0
429         vmask(:,:,:) = 0.0
430      END WHERE
431
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.