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

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

ice sheet coupling: minor changes before merge
File size: 20.7 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   !!
33   !! * Substitutions 
34#  include "domzgr_substitute.h90" 
35#  include "vectopt_loop_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      INTEGER  ::   inum0
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      CHARACTER(20) :: cfile
58      !!----------------------------------------------------------------------
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
120      IF ( lk_vvl ) THEN
121         fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
122         fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
123         fsdept_b(:,:,:) = fsdept_n(:,:,:)
124         fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
125         hu_b (:,:) = hu(:,:)
126         hv_b (:,:) = hv(:,:)
127         hur_b(:,:) = hur(:,:)
128         hvr_b(:,:) = hvr(:,:)
129      END IF
130      !
131   END SUBROUTINE iscpl_stp
132   
133   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
134      !!----------------------------------------------------------------------
135      !!                   ***  ROUTINE iscpl_rst_interpol  ***
136      !!
137      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
138      !!                compute 2d fields of heat, salt and volume correction
139      !!
140      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
141      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
142      !!----------------------------------------------------------------------
143      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
144      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
145      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
146      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
147      !!
148      INTEGER :: ji, jj, jk, iz          !! loop index
149      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
150      !!
151      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
152      REAL(wp):: zdz, zdzm1, zdzp1
153      !!
154      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
155      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
156      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, zhu1, zhv1
157      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
158      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
159      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
160      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
161      !!----------------------------------------------------------------------
162
163      !! allocate variables
164      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
165      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        ) 
166      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
167      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
168      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
169      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
170      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, zhu1, zhv1             ) 
171
172      !! mask value to be sure
173      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
174      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
175     
176      ! compute wmask
177      zwmaskn(:,:,1) = tmask   (:,:,1)
178      zwmaskb(:,:,1) = ptmask_b(:,:,1)
179      DO jk = 2,jpk
180         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
181         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
182      END DO
183           
184      ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
185      sshb (:,:)=sshn(:,:)
186      zssh0(:,:)=sshn(:,:)
187      zsmask0(:,:) = psmask_b(:,:)
188      zsmask1(:,:) = psmask_b(:,:)
189      DO iz = 1,10    ! need to be tuned (configuration dependent) (OK for ISOMIP+)
190         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
191         DO jj = 2,jpj-1
192            DO ji = fs_2, fs_jpim1   ! vector opt.
193               jip1=ji+1; jim1=ji-1;
194               jjp1=jj+1; jjm1=jj-1;
195               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
196               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
197                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
198                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
199                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
200                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
201                  zsmask1(ji,jj)=1._wp
202               END IF
203            END DO
204         END DO
205         CALL lbc_lnk(sshn,'T',1._wp)
206         CALL lbc_lnk(zsmask1,'T',1._wp)
207         zssh0   = sshn
208         zsmask0 = zsmask1
209      END DO
210      sshn(:,:) = sshn(:,:) * ssmask(:,:)
211
212!=============================================================================
213      IF ( lk_vvl ) THEN
214      ! Reconstruction of all vertical scale factors at now time steps
215      ! =============================================================================
216      ! Horizontal scale factor interpolations
217      ! --------------------------------------
218         DO jk = 1,jpk
219            DO jj=1,jpj
220               DO ji=1,jpi
221                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
222                     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) )
223                  ENDIF
224               END DO
225            END DO
226         END DO
227
228         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
229         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
230         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
231
232      ! Vertical scale factor interpolations
233      ! ------------------------------------
234         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
235         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
236         CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
237
238      ! t- and w- points depth
239      ! ----------------------
240         fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
241         fsdepw_n(:,:,1) = 0.0_wp
242         fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
243         DO jj = 1,jpj
244            DO ji = 1,jpi
245               DO jk = 2,mikt(ji,jj)-1
246                  fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
247                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
248                  fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
249               END DO
250               IF (mikt(ji,jj) > 1) THEN
251                  jk = mikt(ji,jj)
252                  fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
253                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
254                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
255               END IF
256               DO jk = mikt(ji,jj)+1, jpk
257                  fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
258                  fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
259                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
260               END DO
261            END DO
262         END DO
263
264      ! t-, u- and v- water column thickness
265      ! ------------------------------------
266         ht(:,:) = 0._wp ; hu(:,:) = 0._wp ; hv(:,:) = 0._wp
267         DO jk = 1, jpkm1
268            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
269            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
270            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
271         END DO
272         !                                        ! Inverse of the local depth
273         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
274         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
275
276      END IF
277
278!=============================================================================
279! compute velocity
280! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
281      ub(:,:,:)=un(:,:,:)
282      vb(:,:,:)=vn(:,:,:)
283      DO jk = 1,jpk
284         DO jj = 1,jpj
285            DO ji = 1,jpi
286               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);
287               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);
288            END DO
289         END DO
290      END DO
291
292! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
293! compute barotropic velocity now and after
294      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
295      zbub(:,:)   = SUM(ztrp,DIM=3)
296      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
297      zbvb(:,:)   = SUM(ztrp,DIM=3)
298      ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:); 
299      zbun(:,:)   = SUM(ztrp,DIM=3)
300      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
301      zbvn(:,:)   = SUM(ztrp,DIM=3)
302
303      ! new water column
304      zhu1=0.0_wp ;
305      zhv1=0.0_wp ;
306      DO jk  = 1,jpk
307        zhu1(:,:) = zhu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
308        zhv1(:,:) = zhv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
309      END DO
310     
311      ! compute correction     
312      zucorr = 0._wp
313      zvcorr = 0._wp
314      DO jj = 1,jpj
315         DO ji = 1,jpi
316            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
317               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
318            END IF
319            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
320               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
321            END IF
322         END DO
323      END DO 
324     
325      ! update velocity
326      DO jk  = 1,jpk
327         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
328         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
329      END DO
330
331!=============================================================================
332      ! compute temp and salt
333      ! compute new tn and sn if we open a new cell
334      tsb (:,:,:,:) = tsn(:,:,:,:)
335      zts0(:,:,:,:) = tsn(:,:,:,:)
336      ztmask1(:,:,:) = ptmask_b(:,:,:)
337      ztmask0(:,:,:) = ptmask_b(:,:,:)
338      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
339          DO jk = 1,jpk-1
340              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
341              DO jj = 2,jpj-1
342                 DO ji = fs_2,fs_jpim1
343                      jip1=ji+1; jim1=ji-1;
344                      jjp1=jj+1; jjm1=jj-1;
345                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
346                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
347                      !! horizontal basic extrapolation
348                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
349                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
350                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
351                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
352                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
353                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
354                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
355                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
356                         ztmask1(ji,jj,jk)=1
357                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
358                      !! vertical extrapolation if horizontal extrapolation failed
359                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
360                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
361                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
362                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
363                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
364                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
365                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
366                            ztmask1(ji,jj,jk)=1._wp
367                         END IF
368                      END IF
369                  END DO
370              END DO
371          END DO
372         
373          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
374          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
375          CALL lbc_lnk(ztmask1,     'T',1._wp)
376
377          ! update
378          zts0(:,:,:,:) = tsn(:,:,:,:)
379          ztmask0 = ztmask1
380
381      END DO
382
383      ! mask new tsn field
384      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
385      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
386
387      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
388      IF ( lk_vvl ) THEN
389         DO jk = 2,jpk-1
390            DO jj = 1,jpj
391               DO ji = 1,jpi
392                  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
393                     !compute weight
394                     zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
395                     zdz   =           fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
396                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  )  - fsdepw_n(ji,jj,jk  ))
397                     IF (zdz .LT. 0._wp) THEN
398                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
399                     END IF
400                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
401                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
402                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/fse3t_n(ji,jj,jk)
403                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
404                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
405                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/fse3t_n(ji,jj,jk)
406                  END IF
407               END DO
408            END DO
409         END DO               
410      END IF
411
412      ! closed pool
413      ! -----------------------------------------------------------------------------------------
414      ! case we open a cell but no neigbour cells available to get an estimate of T and S
415      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
416         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
417         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
418         umask(:,:,:) = 0._wp
419         vmask(:,:,:) = 0._wp
420      END WHERE
421     
422      ! set mbkt and mikt to 1 in thiese location
423      WHERE (SUM(tmask,dim=3) == 0)
424         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
425         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
426      END WHERE
427      ! -------------------------------------------------------------------------------------------
428      ! compute new tn and sn if we close cell
429      ! nothing to do
430      !
431      ! deallocation tmp arrays
432      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
433      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
434      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d              ) 
435      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
436      CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
437      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
438      CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , zhu1 , zhv1          ) 
439
440   END SUBROUTINE iscpl_rst_interpol
441
442END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.