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/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 @ 9630

Last change on this file since 9630 was 9630, checked in by antsia, 6 years ago

Code for divergence correction
File size: 24.0 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, jtop
150      !!
151      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb, zhdivdiff, zfp_ui, zfm_ui, zfp_vi, zfm_vi
152      REAL(wp):: zdz, zdzm1, zdzp1, zvol
153      REAL(wp),DIMENSION(2)::ztrdtsb
154      !!
155      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
156      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
157      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, zhu1, zhv1
158      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
159      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp, zuflux, zvflux
160      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
161      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
162      !!----------------------------------------------------------------------
163
164      !! allocate variables
165      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
166      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        ) 
167      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb, zuflux, zvflux                       ) 
168      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
169      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
170      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
171      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, zhu1, zhv1             ) 
172
173      !! mask value to be sure
174      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
175      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
176     
177      ! compute wmask
178      zwmaskn(:,:,1) = tmask   (:,:,1)
179      zwmaskb(:,:,1) = ptmask_b(:,:,1)
180      DO jk = 2,jpk
181         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
182         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
183      END DO
184           
185      ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
186      sshb (:,:)=sshn(:,:)
187      zssh0(:,:)=sshn(:,:)
188      zsmask0(:,:) = psmask_b(:,:)
189      zsmask1(:,:) = psmask_b(:,:)
190      DO iz = 1,nn_drown    ! need to be tuned (configuration dependent) (OK for ISOMIP+)
191         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
192         DO jj = 2,jpj-1
193            DO ji = fs_2, fs_jpim1   ! vector opt.
194               jip1=ji+1; jim1=ji-1;
195               jjp1=jj+1; jjm1=jj-1;
196               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
197               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
198                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
199                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
200                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
201                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
202                  zsmask1(ji,jj)=1._wp
203               END IF
204            END DO
205         END DO
206         CALL lbc_lnk(sshn,'T',1._wp)
207         CALL lbc_lnk(zsmask1,'T',1._wp)
208         zssh0   = sshn
209         zsmask0 = zsmask1
210      END DO
211      sshn(:,:) = sshn(:,:) * ssmask(:,:)
212
213!=============================================================================
214      IF ( lk_vvl ) THEN
215      ! Reconstruction of all vertical scale factors at now time steps
216      ! =============================================================================
217      ! Horizontal scale factor interpolations
218      ! --------------------------------------
219         DO jk = 1,jpk
220            DO jj=1,jpj
221               DO ji=1,jpi
222                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
223                     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) )
224                  ENDIF
225               END DO
226            END DO
227         END DO
228
229         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
230         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
231         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
232
233      ! Vertical scale factor interpolations
234      ! ------------------------------------
235         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
236         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
237         CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
238
239      ! t- and w- points depth
240      ! ----------------------
241         fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
242         fsdepw_n(:,:,1) = 0.0_wp
243         fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
244         DO jj = 1,jpj
245            DO ji = 1,jpi
246               DO jk = 2,mikt(ji,jj)-1
247                  fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
248                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
249                  fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
250               END DO
251               IF (mikt(ji,jj) > 1) THEN
252                  jk = mikt(ji,jj)
253                  fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
254                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
255                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
256               END IF
257               DO jk = mikt(ji,jj)+1, jpk
258                  fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
259                  fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
260                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
261               END DO
262            END DO
263         END DO
264
265      ! t-, u- and v- water column thickness
266      ! ------------------------------------
267         ht(:,:) = 0._wp ; hu(:,:) = 0._wp ; hv(:,:) = 0._wp
268         DO jk = 1, jpkm1
269            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
270            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
271            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
272         END DO
273         !                                        ! Inverse of the local depth
274         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
275         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
276
277      END IF
278
279!=============================================================================
280! compute velocity
281
282      ub(:,:,:)=un(:,:,:)
283      vb(:,:,:)=vn(:,:,:)
284
285      un(:,:,:)=ub(:,:,:)*umask(:,:,:)
286      vn(:,:,:)=vb(:,:,:)*vmask(:,:,:)
287   
288
289      !-----initialise the divergence and tracer correction------
290      rhdivdiff(:,:,:) = 0._wp
291      rhdivdiff_trac(:,:,:,:) = 0._wp     
292
293      !-----the change of flux in each velocity cell
294      zuflux(:,:,:) = (pe3u_b(:,:,:)*pumask_b(:,:,:) - fse3u_n(:,:,:)*umask(:,:,:))*ub(:,:,:)
295      zvflux(:,:,:) = (pe3v_b(:,:,:)*pvmask_b(:,:,:) - fse3v_n(:,:,:)*vmask(:,:,:))*vb(:,:,:)
296     
297
298      ! compute divergence and tracer correction (whether the horizontal velocity cell is closed or open)   
299      DO jj = 1,jpj
300         DO ji = 1,jpim1   ! jpim1 rather than jpi for the correct lbc_lnk synchronisation
301            !--------u-cell is closed----------- 
302            DO jk  = 1,jpk
303                 !------the contribution to the west side of u-cell
304                 IF(tmask(ji,jj,jk) == 1) THEN   ! the west cell is still active
305                    zvol = e1e2t(ji,jj)*fse3t(ji,jj,jk)
306                    rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + zuflux(ji,jj,jk)*e2u(ji,jj)/zvol 
307                    rhdivdiff_trac(ji,jj,jk,1:2) =  rhdivdiff_trac(ji,jj,jk,1:2) - zuflux(ji,jj,jk)*tsb(ji,jj,jk,1:2)*e2u(ji,jj)/zvol
308                 ELSEIF(mikt(ji,jj) > 1) THEN  ! the west cell becomes dry, put the waterflux to the top of the west column
309                    jtop=mikt(ji,jj)
310                    zvol=e1e2t(ji,jj)*fse3t(ji,jj,jtop)
311                    rhdivdiff(ji,jj,jtop) = rhdivdiff(ji,jj,jtop) + zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
312                    rhdivdiff_trac(ji,jj,jtop,1:2) = rhdivdiff_trac(ji,jj,jtop,1:2)-zuflux(ji,jj,jk)*tsb(ji,jj,jtop,1:2)*e2u(ji,jj)/zvol
313                 !ELSE ! for the closing T-column : nothing to be done
314                 ENDIF
315
316                 !------the contribution to the east side of u-cell
317                 IF(tmask(ji+1,jj,jk) == 1) THEN  ! the east cell is still active
318                    zvol = e1e2t(ji+1,jj)*fse3t(ji+1,jj,jk)
319                    rhdivdiff(ji+1,jj,jk) = rhdivdiff(ji+1,jj,jk) - zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
320                    rhdivdiff_trac(ji+1,jj,jk,1:2) = rhdivdiff_trac(ji+1,jj,jk,1:2) + zuflux(ji,jj,jk)*tsb(ji+1,jj,jk,1:2)*e2u(ji,jj)/zvol
321                 ELSEIF(mikt(ji+1,jj) > 1) THEN    ! the east cell becomes dry, put the waterflux to the top of the east column
322                    jtop = mikt(ji+1,jj)
323                    zvol = e1e2t(ji+1,jj)*fse3t(ji+1,jj,jtop)
324                    rhdivdiff(ji+1,jj,jtop) = rhdivdiff(ji+1,jj,jtop) - zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
325                    rhdivdiff_trac(ji+1,jj,jtop,1:2) = rhdivdiff_trac(ji+1,jj,jtop,1:2) + zuflux(ji,jj,jk)*tsb(ji+1,jj,jtop,1:2)*e2u(ji,jj)/zvol
326                 !ELSE  ! for the closing T-column : nothing to be done
327                 ENDIF
328            ENDDO
329         END DO
330      END DO
331
332
333      DO jj = 1,jpjm1  ! jpim1 rather than jpi for the correct lbc_lnk synchronisation
334         DO ji = 1,jpi   
335            !-----------v-cell is closed-------------
336            DO jk  = 1,jpk
337                 !------the contribution to the south side of v-cell
338                 IF(tmask(ji,jj,jk) == 1) THEN  ! the south cell is still active
339                   zvol = e1e2t(ji,jj)*fse3t(ji,jj,jk)
340                   rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
341                   rhdivdiff_trac(ji,jj,jk,1:2) =  rhdivdiff_trac(ji,jj,jk,1:2) - zvflux(ji,jj,jk)*tsb(ji,jj,jk,1:2)*e1v(ji,jj)/zvol
342                 ELSEIF(mikt(ji,jj) > 1) THEN  !the south cell becomes dry, put the waterflux to the top of the south column
343                   jtop = mikt(ji,jj)
344                   zvol = e1e2t(ji,jj)*fse3t(ji,jj,jtop)
345                   rhdivdiff(ji,jj,jtop) = rhdivdiff(ji,jj,jtop) + zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
346                   rhdivdiff_trac(ji,jj,jtop,1:2) = rhdivdiff_trac(ji,jj,jtop,1:2) - zvflux(ji,jj,jk)*tsb(ji,jj,jtop,1:2)*e1v(ji,jj)/zvol
347                 ENDIF
348
349                 !------the contribution to the north side of v-cell
350                 IF(tmask(ji,jj+1,jk) == 1) THEN ! the north cell is still active
351                    zvol = e1e2t(ji,jj+1)*fse3t(ji,jj+1,jk)
352                    rhdivdiff(ji,jj+1,jk) = rhdivdiff(ji,jj+1,jk) - zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
353                    rhdivdiff_trac(ji,jj+1,jk,1:2) = rhdivdiff_trac(ji,jj+1,jk,1:2) + zvflux(ji,jj,jk)*tsb(ji,jj+1,jk,1:2)*e1v(ji,jj)/zvol
354                 ELSEIF(mikt(ji,jj+1) > 1) THEN  ! the north cell becomes dry, put the waterflux to the top of the north column
355                    jtop = mikt(ji,jj+1)
356                    zvol = e1e2t(ji,jj+1)*fse3t(ji,jj+1,jtop)
357                    rhdivdiff(ji,jj+1,jtop) = rhdivdiff(ji,jj+1,jtop) - zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
358                    rhdivdiff_trac(ji,jj+1,jtop,1:2) = rhdivdiff_trac(ji,jj+1,jtop,1:2) + zvflux(ji,jj,jk)*tsb(ji,jj+1,jtop,1:2)*e1v(ji,jj)/zvol
359                 ENDIF
360            ENDDO
361
362         END DO
363      END DO
364     
365
366      CALL lbc_lnk(rhdivdiff,   'T',1._wp)
367      CALL lbc_lnk(rhdivdiff_trac(:,:,:,1),   'T',1._wp)
368      CALL lbc_lnk(rhdivdiff_trac(:,:,:,2),   'T',1._wp)
369
370
371      ! compute temp and salt
372      ! compute new tn and sn if we open a new cell
373      tsb (:,:,:,:) = tsn(:,:,:,:)
374      zts0(:,:,:,:) = tsn(:,:,:,:)
375      ztmask1(:,:,:) = ptmask_b(:,:,:)
376      ztmask0(:,:,:) = ptmask_b(:,:,:)
377      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
378          DO jk = 1,jpk-1
379              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
380              DO jj = 2,jpj-1
381                 DO ji = fs_2,fs_jpim1
382                      jip1=ji+1; jim1=ji-1;
383                      jjp1=jj+1; jjm1=jj-1;
384                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
385                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
386                      !! horizontal basic extrapolation
387                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
388                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
389                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
390                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
391                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
392                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
393                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
394                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
395                         ztmask1(ji,jj,jk)=1
396                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
397                      !! vertical extrapolation if horizontal extrapolation failed
398                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
399                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
400                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
401                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
402                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
403                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
404                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
405                            ztmask1(ji,jj,jk)=1._wp
406                         END IF
407                      END IF
408                  END DO
409              END DO
410          END DO
411         
412          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
413          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
414          CALL lbc_lnk(ztmask1,     'T',1._wp)
415
416          ! update
417          zts0(:,:,:,:) = tsn(:,:,:,:)
418          ztmask0 = ztmask1
419
420      END DO
421
422      ! mask new tsn field
423      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
424      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
425
426      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
427!      IF ( lk_vvl ) THEN
428!         DO jk = 2,jpk-1
429!            DO jj = 1,jpj
430!               DO ji = 1,jpi
431!                  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
432!                     !compute weight
433!                     zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw_n(ji,jj,jk+1))
434!                     zdzm1 = MAX(0._wp,gdepw_n(ji,jj,jk  ) - pdepw_b(ji,jj,jk  ))
435!                     zdz   = e3t_n(ji,jj,jk) - zdzp1 - zdzm1 ! if isf : e3t = gdepw_n(ji,jj,jk+1)- gdepw_n(ji,jj,jk)
436!                     !
437!                     IF (zdz .LT. 0._wp) THEN
438!                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
439!                     END IF
440!                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem) &
441!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem) &
442!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem) )/e3t_n(ji,jj,jk)
443!                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal) &
444!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal) &
445!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal) )/e3t_n(ji,jj,jk)
446!                  END IF
447!               END DO
448!            END DO
449!         END DO               
450!      END IF
451
452      ! closed pool
453      ! -----------------------------------------------------------------------------------------
454      ! case we open a cell but no neigbour cells available to get an estimate of T and S
455      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
456         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
457         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
458         umask(:,:,:) = 0._wp
459         vmask(:,:,:) = 0._wp
460      END WHERE
461     
462      ! set mbkt and mikt to 1 in thiese location
463      WHERE (SUM(tmask,dim=3) == 0)
464         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
465         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
466      END WHERE
467
468      ! -------------------------------------------------------------------------------------------
469      ! compute new tn and sn if we close cell
470      ! nothing to do
471      !
472      ! deallocation tmp arrays
473      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
474      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
475      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d, zuflux, zvflux    ) 
476      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
477      CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
478      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
479      CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , zhu1 , zhv1          ) 
480
481   END SUBROUTINE iscpl_rst_interpol
482
483END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.