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.
isfcpl.F90 in NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/ISF – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/ISF/isfcpl.F90 @ 12342

Last change on this file since 12342 was 12342, checked in by acc, 4 years ago

Branch 2019/dev_r12072_MERGE_OPTION2_2019. Implement recommended changes to fix ticket #2370
File size: 34.8 KB
Line 
1MODULE isfcpl
2   !!======================================================================
3   !!                       ***  MODULE  isfcpl  ***
4   !!
5   !! iceshelf coupling module : module managing the coupling between NEMO and an ice sheet model
6   !!
7   !!======================================================================
8   !! History :  4.1  !  2019-07  (P. Mathiot) Original code
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   isfrst : read/write iceshelf variables in/from restart
13   !!----------------------------------------------------------------------
14   USE isf_oce                          ! ice shelf variable
15   USE isfutils, ONLY : debug
16   USE lib_mpp , ONLY: mpp_sum, mpp_max ! mpp routine
17   USE domvvl  , ONLY: dom_vvl_zgr      ! vertical scale factor interpolation
18   USE domngb  , ONLY: dom_ngb          ! find the closest grid point from a given lon/lat position
19   !
20   USE oce            ! ocean dynamics and tracers
21   USE in_out_manager ! I/O manager
22   USE iom            ! I/O library
23   !
24   IMPLICIT NONE
25
26   PRIVATE
27
28   PUBLIC isfcpl_rst_write, isfcpl_init                    ! iceshelf restart read and write
29   PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons  ! iceshelf correction for ssh, tra, dyn and conservation
30
31   TYPE isfcons
32      INTEGER :: ii     ! i global
33      INTEGER :: jj     ! j global
34      INTEGER :: kk     ! k level
35      REAL(wp):: dvol   ! volume increment
36      REAL(wp):: dsal   ! salt increment
37      REAL(wp):: dtem   ! heat increment
38      REAL(wp):: lon    ! lon
39      REAL(wp):: lat    ! lat
40      INTEGER :: ngb    ! 0/1 (valid location or not (ie on halo or no neigbourg))
41   END TYPE
42   !
43   !!----------------------------------------------------------------------
44   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
45   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
46   !! Software governed by the CeCILL license (see ./LICENSE)
47   !!----------------------------------------------------------------------
48CONTAINS
49   SUBROUTINE isfcpl_init()
50      !!---------------------------------------------------------------------
51      !!                   ***  ROUTINE iscpl_init  ***
52      !!
53      !! ** Purpose : correct ocean state for new wet cell and horizontal divergence
54      !!              correction for the dynamical adjustement
55      !!
56      !! ** Action : - compute ssh on new wet cell
57      !!             - compute T/S on new wet cell
58      !!             - compute horizontal divergence correction as a volume flux
59      !!             - compute the T/S/vol correction increment to keep trend to 0
60      !!
61      !!---------------------------------------------------------------------
62      INTEGER :: id
63      !!----------------------------------------------------------------------
64      !
65      ! start on an euler time step
66      neuler = 0
67      !
68      ! allocation and initialisation to 0
69      CALL isf_alloc_cpl()
70      !
71      ! check presence of variable needed for coupling
72      ! iom_varid return 0 if not found
73      id = 1
74      id = id * iom_varid(numror, 'ssmask', ldstop = .false.)
75      id = id * iom_varid(numror, 'tmask' , ldstop = .false.)
76      id = id * iom_varid(numror, 'e3t_n' , ldstop = .false.)
77      id = id * iom_varid(numror, 'e3u_n' , ldstop = .false.)
78      id = id * iom_varid(numror, 'e3v_n' , ldstop = .false.)
79      IF(lwp) WRITE(numout,*) ' isfcpl_init:', id
80      IF (id == 0) THEN
81         IF(lwp) WRITE(numout,*) ' isfcpl_init: restart variables for ice sheet coupling are missing, skip coupling for this leg ' 
82         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
83         IF(lwp) WRITE(numout,*) ''
84      ELSE
85         ! extrapolation ssh
86         CALL isfcpl_ssh()
87         !
88         ! extrapolation tracer properties
89         CALL isfcpl_tra()
90         !
91         ! correction of the horizontal divergence and associated temp. and salt content flux
92         ! Need to : - include in the cpl cons the risfcpl_vol/tsc contribution
93         !           - decide how to manage thickness level change in conservation
94         CALL isfcpl_vol()
95         !
96         ! apply the 'conservation' method
97         IF ( ln_isfcpl_cons ) CALL isfcpl_cons()
98         !
99      END IF
100      !
101      ! mask velocity properly (mask used in restart not compatible with new mask)
102      un(:,:,:) = un(:,:,:) * umask(:,:,:)
103      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
104      !
105      ! all before fields set to now values
106      tsb  (:,:,:,:) = tsn  (:,:,:,:)
107      ub   (:,:,:)   = un   (:,:,:)
108      vb   (:,:,:)   = vn   (:,:,:)
109      sshb (:,:)     = sshn (:,:)
110      e3t_b(:,:,:)   = e3t_n(:,:,:)
111 
112      ! prepare writing restart
113      IF( lwxios ) THEN
114         CALL iom_set_rstw_var_active('ssmask')
115         CALL iom_set_rstw_var_active('tmask')
116         CALL iom_set_rstw_var_active('e3t_n')
117         CALL iom_set_rstw_var_active('e3u_n')
118         CALL iom_set_rstw_var_active('e3v_n')
119      END IF
120      !
121   END SUBROUTINE isfcpl_init
122   !
123   SUBROUTINE isfcpl_rst_write(kt)
124      !!---------------------------------------------------------------------
125      !!                   ***  ROUTINE iscpl_rst_write  ***
126      !!
127      !! ** Purpose : write icesheet coupling variables in restart
128      !!
129      !!-------------------------- IN  --------------------------------------
130      INTEGER, INTENT(in) :: kt
131      !!----------------------------------------------------------------------
132      !
133      IF( lwxios ) CALL iom_swap( cwxios_context )
134      CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask , ldxios = lwxios )
135      CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask, ldxios = lwxios )
136      CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , e3t_n , ldxios = lwxios )
137      CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , e3u_n , ldxios = lwxios )
138      CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , e3v_n , ldxios = lwxios )
139      CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n , ldxios = lwxios )
140      IF( lwxios ) CALL iom_swap( cxios_context )
141      !
142   END SUBROUTINE isfcpl_rst_write
143
144   SUBROUTINE isfcpl_ssh()
145      !!----------------------------------------------------------------------
146      !!                   ***  ROUTINE iscpl_ssh  ***
147      !!
148      !! ** Purpose :   basic guess of ssh in new wet cell
149      !!
150      !! ** Method  :   basic extrapolation from neigbourg cells
151      !!
152      !!----------------------------------------------------------------------
153      !!
154      INTEGER :: ji, jj, jd, jk      !! loop index
155      INTEGER :: jip1, jim1, jjp1, jjm1
156      !!
157      REAL(wp):: zsummsk
158      REAL(wp), DIMENSION(jpi,jpj) :: zdssmask, zssmask0, zssmask_b, zssh
159      !!----------------------------------------------------------------------
160      !
161      CALL iom_get( numror, jpdom_autoglo, 'ssmask'  , zssmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
162
163      ! compute new ssh if we open a full water column
164      ! rude average of the closest neigbourgs (e1e2t not taking into account)
165      !
166      zssh(:,:)     = sshn(:,:)
167      zssmask0(:,:) = zssmask_b(:,:)
168      !
169      DO jd = 1, nn_drown
170         !
171         zdssmask(:,:) = ssmask(:,:) - zssmask0(:,:)
172         DO jj = 2,jpj-1
173            DO ji = 2,jpi-1
174               jip1=ji+1; jim1=ji-1;
175               jjp1=jj+1; jjm1=jj-1;
176               !
177               zsummsk = zssmask0(jip1,jj) + zssmask0(jim1,jj) + zssmask0(ji,jjp1) + zssmask0(ji,jjm1)
178               !
179               IF (zdssmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp) THEN
180                  sshn(ji,jj)=( zssh(jip1,jj)*zssmask0(jip1,jj)     &
181                  &           + zssh(jim1,jj)*zssmask0(jim1,jj)     &
182                  &           + zssh(ji,jjp1)*zssmask0(ji,jjp1)     &
183                  &           + zssh(ji,jjm1)*zssmask0(ji,jjm1)) / zsummsk
184                  zssmask_b(ji,jj) = 1._wp
185               ENDIF
186            END DO
187         END DO
188         !
189         zssh(:,:) = sshn(:,:)
190         zssmask0(:,:) = zssmask_b(:,:)
191         !
192         CALL lbc_lnk_multi( 'iscplrst', zssh, 'T', 1., zssmask0, 'T', 1. )
193         !
194      END DO
195      !
196      ! update sshn
197      sshn(:,:) = zssh(:,:) * ssmask(:,:)
198      !
199      sshb(:,:) = sshn(:,:)
200      !
201      IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',sshn(:,:))
202      !
203      ! recompute the vertical scale factor, depth and water thickness
204      IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from sshn (new wet cell)'
205      IF(lwp) write(numout,*) '~~~~~~~~~~~'
206      DO jk = 1, jpk
207         e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
208             &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
209             &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
210      END DO
211      e3t_b(:,:,:) = e3t_n(:,:,:)
212      CALL dom_vvl_zgr()
213      !
214   END SUBROUTINE isfcpl_ssh
215
216   SUBROUTINE isfcpl_tra()
217      !!----------------------------------------------------------------------
218      !!                   ***  ROUTINE iscpl_tra  ***
219      !!
220      !! ** Purpose :   compute new tn, sn in case of evolving geometry of ice shelves
221      !!
222      !! ** Method  :   tn, sn : basic extrapolation from neigbourg cells
223      !!
224      !!----------------------------------------------------------------------
225      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
226      !REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
227      !!
228      INTEGER :: ji, jj, jk, jd          !! loop index
229      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
230      !!
231      REAL(wp):: zsummsk
232      REAL(wp):: zdz, zdzm1, zdzp1
233      !!
234      REAL(wp), DIMENSION(jpi,jpj)          :: zdmask 
235      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn
236      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
237      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
238      !!----------------------------------------------------------------------
239      !
240      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
241      !CALL iom_get( numror, jpdom_autoglo, 'wmask'  , zwmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
242      !CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
243      !
244      !
245      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
246      !PM: Is this IF needed since change to VVL by default
247      !bugged : to be corrected (PM)
248      ! back up original t/s/mask
249      !tsb (:,:,:,:) = tsn(:,:,:,:)
250      !
251     ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
252
253!      IF (.NOT.ln_linssh) THEN
254!         DO jk = 2,jpk-1
255!            DO jj = 1,jpj
256!               DO ji = 1,jpi
257!                  IF (wmask(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ztmask_b(ji,jj,1)==0._wp) ) THEN
258!
259!                     !compute weight
260!                     zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw_n(ji,jj,jk+1))
261!                     zdzm1 = MAX(0._wp,gdepw_n(ji,jj,jk  ) - pdepw_b(ji,jj,jk  ))
262!                     zdz   = e3t_n(ji,jj,jk) - zdzp1 - zdzm1 ! if isf : e3t = gdepw_n(ji,jj,jk+1)- gdepw_n(ji,jj,jk)
263!
264!                     IF (zdz .LT. 0._wp) THEN
265!                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
266!                     END IF
267!
268!                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem) &
269!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem) &
270!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem) )/e3t_n(ji,jj,jk)
271!
272!                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal) &
273!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal) &
274!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal) )/e3t_n(ji,jj,jk)
275!
276!                  END IF
277!               END DO
278!            END DO
279!         END DO
280!      END IF
281
282      zts0(:,:,:,:)  = tsn(:,:,:,:)
283      ztmask0(:,:,:) = ztmask_b(:,:,:)
284      ztmask1(:,:,:) = ztmask_b(:,:,:)
285      !
286      ! iterate the extrapolation processes nn_drown times
287      DO jd = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
288         DO jk = 1,jpk-1
289            !
290            ! define new wet cell
291            zdmask(:,:) = tmask(:,:,jk) - ztmask0(:,:,jk);
292            !
293            DO jj = 2,jpj-1
294               DO ji = 2,jpi-1
295                  jip1=ji+1; jim1=ji-1;
296                  jjp1=jj+1; jjm1=jj-1;
297                  !
298                  ! check if a wet neigbourg cell is present
299                  zsummsk = ztmask0(jip1,jj  ,jk) + ztmask0(jim1,jj  ,jk) &
300                          + ztmask0(ji  ,jjp1,jk) + ztmask0(ji  ,jjm1,jk)
301                  !
302                  ! if neigbourg wet cell available at the same level
303                  IF ( zdmask(ji,jj) == 1._wp  .AND. zsummsk /= 0._wp ) THEN
304                     !
305                     ! horizontal basic extrapolation
306                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1) * ztmask0(jip1,jj  ,jk) &
307                     &               + zts0(jim1,jj  ,jk,1) * ztmask0(jim1,jj  ,jk) &
308                     &               + zts0(ji  ,jjp1,jk,1) * ztmask0(ji  ,jjp1,jk) &
309                     &               + zts0(ji  ,jjm1,jk,1) * ztmask0(ji  ,jjm1,jk) ) / zsummsk
310                     tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2) * ztmask0(jip1,jj  ,jk) &
311                     &               + zts0(jim1,jj  ,jk,2) * ztmask0(jim1,jj  ,jk) &
312                     &               + zts0(ji  ,jjp1,jk,2) * ztmask0(ji  ,jjp1,jk) &
313                     &               + zts0(ji  ,jjm1,jk,2) * ztmask0(ji  ,jjm1,jk) ) / zsummsk
314                     !
315                     ! update mask for next pass
316                     ztmask1(ji,jj,jk)=1
317                     !
318                  ! in case no neigbourg wet cell available at the same level
319                  ! check if a wet cell is available below
320                  ELSEIF (zdmask(ji,jj) == 1._wp .AND. zsummsk == 0._wp) THEN
321                     !
322                     ! vertical extrapolation if horizontal extrapolation failed
323                     jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
324                     !
325                     ! check if a wet neigbourg cell is present
326                     zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
327                     IF (zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
328                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
329                        &               + zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / zsummsk
330                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
331                        &               + zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / zsummsk
332                        !
333                        ! update mask for next pass
334                        ztmask1(ji,jj,jk)=1._wp
335                     END IF
336                  END IF
337               END DO
338            END DO
339         END DO
340         !
341         ! update temperature and salinity and mask
342         zts0(:,:,:,:)  = tsn(:,:,:,:)
343         ztmask0(:,:,:) = ztmask1(:,:,:)
344         !
345         CALL lbc_lnk_multi( 'iscplrst', zts0(:,:,:,jp_tem), 'T', 1., zts0(:,:,:,jp_sal), 'T', 1., ztmask0, 'T', 1.)
346         !
347      END DO  ! nn_drown
348      !
349      ! mask new tsn field
350      tsn(:,:,:,jp_tem) = zts0(:,:,:,jp_tem) * tmask(:,:,:)
351      tsn(:,:,:,jp_sal) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
352      !
353      ! sanity check
354      ! -----------------------------------------------------------------------------------------
355      ! case we open a cell but no neigbour cells available to get an estimate of T and S
356      DO jk = 1,jpk-1
357         DO jj = 1,jpj
358            DO ji = 1,jpi
359               IF (tmask(ji,jj,jk) == 1._wp .AND. tsn(ji,jj,jk,2) == 0._wp)              &
360                  &   CALL ctl_stop('STOP', 'failing to fill all new weet cell,     &
361                  &                          try increase nn_drown or activate XXXX &
362                  &                         in your domain cfg computation'         )
363            END DO
364         END DO
365      END DO
366      !
367   END SUBROUTINE isfcpl_tra
368
369   SUBROUTINE isfcpl_vol()
370      !!----------------------------------------------------------------------
371      !!                   ***  ROUTINE iscpl_vol  ***
372      !!
373      !! ** Purpose : compute the correction of the local divergence to apply 
374      !!              during the first time step after the coupling.
375      !!
376      !! ** Method  : - compute horizontal vol div. before/after coupling
377      !!              - compute vertical input
378      !!              - compute correction
379      !!               
380      !!----------------------------------------------------------------------
381      !!
382      INTEGER :: ji, jj, jk 
383      INTEGER :: ikb, ikt
384      !!
385      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zqvolb, zqvoln  ! vol flux div.         before/after coupling
386      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3u_b, ze3v_b  ! vertical scale factor before/after coupling
387      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b        ! mask                  before       coupling
388      !!----------------------------------------------------------------------
389      !
390      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios )
391      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b  , ldxios = lrxios )
392      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b  , ldxios = lrxios )
393      !
394      ! 1.0: compute horizontal volume flux divergence difference before-after coupling
395      !
396      DO jk = 1, jpk                                 ! Horizontal slab
397         ! 1.1: get volume flux before coupling (>0 out)
398         DO jj = 2, jpjm1
399            DO ji = 2, jpim1
400               zqvolb(ji,jj,jk) =  (   e2u(ji,jj) * ze3u_b(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  ) * ze3u_b(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
401                  &                  + e1v(ji,jj) * ze3v_b(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * ze3v_b(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  ) &
402                  &                * ztmask_b(ji,jj,jk)
403            END DO
404         ENDDO
405         !
406         ! 1.2: get volume flux after coupling (>0 out)
407         ! properly mask velocity
408         ! (velocity are still mask with old mask at this stage)
409         un(:,:,jk) = un(:,:,jk) * umask(:,:,jk)
410         vn(:,:,jk) = vn(:,:,jk) * vmask(:,:,jk)
411         ! compute volume flux divergence after coupling
412         DO jj = 2, jpjm1
413            DO ji = 2, jpim1
414               zqvoln(ji,jj,jk) = (   e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  ) * e3u_n(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
415                  &                 + e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * e3v_n(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  ) &
416                  &               * tmask(ji,jj,jk)
417            END DO
418         ENDDO
419         !
420         ! 1.3: get 3d volume flux difference (before - after cpl) (>0 out)
421         !      correction to add is _b - _n
422         risfcpl_vol(:,:,jk) = zqvolb(:,:,jk) - zqvoln(:,:,jk)
423      END DO
424      !
425      ! 2.0: include the contribution of the vertical velocity in the volume flux correction
426      !
427      DO jj = 2, jpjm1
428         DO ji = 2, jpim1
429            !
430            ikt = mikt(ji,jj)
431            IF ( ikt > 1 .AND. ssmask(ji,jj) == 1 ) THEN
432               risfcpl_vol(ji,jj,ikt) = risfcpl_vol(ji,jj,ikt) + SUM(zqvolb(ji,jj,1:ikt-1))  ! test sign
433            ENDIF
434            !
435         END DO
436      ENDDO
437      !
438      CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. )
439      !
440      ! 3.0: set total correction (div, tra, ssh)
441      !
442      ! 3.1: mask volume flux divergence correction
443      risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
444      !
445      ! 3.2: get 3d tra increment to apply at the first time step
446      ! temperature and salt content flux computed using local tsn
447      ! (very simple advection scheme)
448      ! (>0 out)
449      risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_tem)
450      risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_sal)
451      !
452      ! 3.3: ssh correction (for dynspg_ts)
453      risfcpl_ssh(:,:) = 0.0
454      DO jk = 1,jpk
455         risfcpl_ssh(:,:) = risfcpl_ssh(:,:) + risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
456      END DO
457
458   END SUBROUTINE isfcpl_vol
459
460   SUBROUTINE isfcpl_cons()
461      !!----------------------------------------------------------------------
462      !!                   ***  ROUTINE iscpl_cons  ***
463      !!
464      !! ** Purpose :   compute the corrective increment in volume/salt/heat to put back the vol/heat/salt
465      !!                removed or added during the coupling processes (wet or dry new cell)
466      !!
467      !! ** Method  :   - compare volume/heat/salt before and after
468      !!                - look for the closest wet cells (share amoung neigbourgs if there are)
469      !!                - build the correction increment to applied at each time step
470      !!               
471      !!----------------------------------------------------------------------
472      !
473      TYPE(isfcons), DIMENSION(:),ALLOCATABLE :: zisfpts ! list of point receiving a correction
474      !
475      INTEGER  ::   ji   , jj  , jk  , jproc          ! loop index
476      INTEGER  ::   jip1 , jim1, jjp1, jjm1           ! dummy indices
477      INTEGER  ::   iig  , ijg, ik                    ! dummy indices
478      INTEGER  ::   jisf                              ! start, end and current position in the increment array
479      INTEGER  ::   ingb, ifind, nisfl_area           ! global number of cell concerned by the wet->dry case
480      INTEGER, DIMENSION(jpnij) :: nisfl              ! local  number of cell concerned by the wet->dry case
481      !
482      REAL(wp) ::   z1_sum, z1_rdtiscpl
483      REAL(wp) ::   zdtem, zdsal, zdvol, zratio       ! tem, sal, vol increment
484      REAL(wp) ::   zlon , zlat                       ! target location 
485      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b    ! mask before
486      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b      ! scale factor before
487      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zt_b      ! scale factor before
488      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zs_b      ! scale factor before
489      !!----------------------------------------------------------------------
490
491      !==============================================================================
492      ! 1.0: initialisation
493      !==============================================================================
494
495      ! get restart variable
496      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b(:,:,:), ldxios = lrxios   ) ! need to extrapolate T/S
497      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:)  , ldxios = lrxios )
498      CALL iom_get( numror, jpdom_autoglo, 'tn'     , zt_b(:,:,:)    , ldxios = lrxios )
499      CALL iom_get( numror, jpdom_autoglo, 'sn'     , zs_b(:,:,:)    , ldxios = lrxios )
500
501      ! compute run length
502      nstp_iscpl  = nitend - nit000 + 1
503      rdt_iscpl   = nstp_iscpl * rn_rdt
504      z1_rdtiscpl = 1._wp / rdt_iscpl 
505
506      IF (lwp) WRITE(numout,*) '            nb of stp for cons  = ', nstp_iscpl
507      IF (lwp) WRITE(numout,*) '            coupling time step  = ', rdt_iscpl
508
509      ! initialisation correction
510      risfcpl_cons_vol = 0.0
511      risfcpl_cons_ssh = 0.0
512      risfcpl_cons_tsc = 0.0
513
514      !==============================================================================
515      ! 2.0: diagnose the heat, salt and volume input and compute the correction variable
516      !      for case where we wet a cell or cell still wet (no change in cell status)
517      !==============================================================================
518
519      DO jk = 1,jpk-1
520         DO jj = nldj,nlej
521            DO ji = nldi,nlei
522
523               ! volume diff
524               zdvol = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
525
526               ! heat diff
527               zdtem = tsn (ji,jj,jk,jp_tem) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
528                     - zt_b(ji,jj,jk)        * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
529
530               ! salt diff
531               zdsal = tsn(ji,jj,jk,jp_sal) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
532                     - zs_b(ji,jj,jk)       * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
533           
534               ! volume, heat and salt differences in each cell (>0 means correction is an outward flux)
535               ! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
536               risfcpl_cons_vol(ji,jj,jk)        = (   zdvol * e1e2t(ji,jj) + risfcpl_vol(ji,jj,jk)        ) * z1_rdtiscpl
537               risfcpl_cons_tsc(ji,jj,jk,jp_sal) = ( - zdsal * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_sal) ) * z1_rdtiscpl
538               risfcpl_cons_tsc(ji,jj,jk,jp_tem) = ( - zdtem * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_tem) ) * z1_rdtiscpl
539
540            END DO
541         END DO
542      END DO
543      !
544      !==============================================================================
545      ! 3.0: diagnose the heat, salt and volume input and compute the correction variable
546      !      for case where we close a cell
547      !==============================================================================
548      !
549      ! compute the total number of point receiving a correction increment for each processor
550      ! local
551      nisfl(:)=0
552      DO jk = 1,jpk-1
553         DO jj = nldj,nlej
554            DO ji = nldi,nlei
555               jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
556               IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
557            ENDDO
558         ENDDO
559      ENDDO
560      !
561      ! global
562      CALL mpp_sum('isfcpl',nisfl  )
563      !
564      ! allocate list of point receiving correction
565      ALLOCATE(zisfpts(nisfl(narea)))
566      !
567      zisfpts(:) = isfcons(0,0,0,-HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), 0)
568      !
569      ! start computing the correction and fill zisfpts
570      ! local
571      jisf = 0
572      DO jk = 1,jpk-1
573         DO jj = nldj,nlej
574            DO ji = nldi,nlei
575               IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
576
577                  jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
578
579                  zdvol = risfcpl_cons_vol(ji,jj,jk       )
580                  zdsal = risfcpl_cons_tsc(ji,jj,jk,jp_sal)
581                  zdtem = risfcpl_cons_tsc(ji,jj,jk,jp_tem)
582
583                  IF ( SUM( tmask(jim1:jip1,jjm1:jjp1,jk) ) > 0._wp ) THEN
584                     ! spread correction amoung neigbourg wet cells (horizontal direction first)
585                     ! as it is a rude correction corner and lateral cell have the same weight
586                     !
587                     z1_sum =  1._wp / SUM( tmask(jim1:jip1,jjm1:jjp1,jk) )
588                     !
589                     ! lateral cells
590                     IF (tmask(jip1,jj  ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jj  , jk, zdvol, zdsal, zdtem, z1_sum)
591                     IF (tmask(jim1,jj  ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jj  , jk, zdvol, zdsal, zdtem, z1_sum)
592                     IF (tmask(ji  ,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji  , jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
593                     IF (tmask(ji  ,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji  , jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
594                     !
595                     ! corner  cells
596                     IF (tmask(jip1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
597                     IF (tmask(jim1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
598                     IF (tmask(jim1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
599                     IF (tmask(jip1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
600                     !
601                  ELSE IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
602                     ! spread correction amoung neigbourg wet cells (vertical direction)
603                     CALL update_isfpts(zisfpts, jisf, ji  , jj  , jk+1, zdvol, zdsal, zdtem, 1., 0)
604                  ELSE
605                     ! need to find where to put correction in later on
606                     CALL update_isfpts(zisfpts, jisf, ji  , jj  , jk  , zdvol, zdsal, zdtem, 1., 1)
607                  END IF
608               END IF
609            END DO
610         END DO
611      END DO
612      !
613      ! share data among all processes because for some point we need to find the closest wet point (could be on other process)
614      DO jproc=1,jpnij
615         !
616         ! share total number of isf point treated for proc jproc
617         IF (jproc==narea) THEN
618            nisfl_area=nisfl(jproc)
619         ELSE
620            nisfl_area=0
621         END IF
622         CALL mpp_max('isfcpl',nisfl_area)
623         !
624         DO jisf = 1,nisfl_area
625            !
626            IF (jproc==narea) THEN
627               ! indices (conversion to global indices and sharing)
628               iig = zisfpts(jisf)%ii       ; ijg = zisfpts(jisf)%jj       ; ik = zisfpts(jisf)%kk
629               !
630               ! data
631               zdvol = zisfpts(jisf)%dvol   ; zdsal = zisfpts(jisf)%dsal   ; zdtem = zisfpts(jisf)%dtem
632               !
633               ! location
634               zlat = zisfpts(jisf)%lat     ; zlon = zisfpts(jisf)%lon
635               !
636               ! find flag
637               ingb = zisfpts(jisf)%ngb
638            ELSE
639               iig  =0   ; ijg  =0   ; ik   =0 
640               zdvol=-HUGE(1.0) ; zdsal=-HUGE(1.0) ; zdtem=-HUGE(1.0)
641               zlat =-HUGE(1.0) ; zlon =-HUGE(1.0)   
642               ingb = 0
643            END IF
644            !
645            ! share data (need synchronisation of data as get_correction call a global com)
646            CALL mpp_max('isfcpl',iig)   ; CALL mpp_max('isfcpl',ijg)   ; CALL mpp_max('isfcpl',ik)
647            CALL mpp_max('isfcpl',zdvol) ; CALL mpp_max('isfcpl',zdsal) ; CALL mpp_max('isfcpl',zdtem)
648            CALL mpp_max('isfcpl',zlat)  ; CALL mpp_max('isfcpl',zlon)
649            CALL mpp_max('isfcpl',ingb)
650            !
651            ! fill the 3d correction array
652            CALL get_correction(iig, ijg, ik, zlon, zlat, zdvol, zdsal, zdtem, ingb)
653         END DO
654      END DO
655      !
656      !==============================================================================
657      ! 4.0: finalisation and compute ssh equivalent of the volume correction
658      !==============================================================================
659      !
660      ! mask (>0 out)
661      risfcpl_cons_vol(:,:,:       ) = risfcpl_cons_vol(:,:,:       ) * tmask(:,:,:)
662      risfcpl_cons_tsc(:,:,:,jp_sal) = risfcpl_cons_tsc(:,:,:,jp_sal) * tmask(:,:,:)
663      risfcpl_cons_tsc(:,:,:,jp_tem) = risfcpl_cons_tsc(:,:,:,jp_tem) * tmask(:,:,:)
664      !
665      ! add lbclnk
666      CALL lbc_lnk_multi( 'iscplrst', risfcpl_cons_tsc(:,:,:,jp_tem), 'T', 1., risfcpl_cons_tsc(:,:,:,jp_sal), 'T', 1., &
667         &                            risfcpl_cons_vol(:,:,:)       , 'T', 1.)
668      !
669      ! ssh correction (for dynspg_ts)
670      DO jk = 1,jpk
671         risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) + risfcpl_cons_vol(:,:,jk)
672      END DO
673      risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) * r1_e1e2t(:,:)
674      !
675   END SUBROUTINE isfcpl_cons
676   !
677   SUBROUTINE update_isfpts(sisfpts, kpts, ki, kj, kk, pdvol, pdsal, pdtem, pratio, kfind)
678      !!---------------------------------------------------------------------
679      !!                  ***  ROUTINE update_isfpts  ***
680      !!
681      !! ** Purpose : if a cell become dry, we need to put the corrective increment elsewhere
682      !!
683      !! ** Action  : update the list of point
684      !!
685      !!----------------------------------------------------------------------
686      !!----------------------------------------------------------------------
687      TYPE(isfcons), DIMENSION(:), INTENT(inout) :: sisfpts
688      INTEGER,                     INTENT(inout) :: kpts
689      !!----------------------------------------------------------------------
690      INTEGER,      INTENT(in   )           :: ki, kj, kk                  !    target location (kfind=0)
691      !                                                                    ! or source location (kfind=1)
692      INTEGER,      INTENT(in   ), OPTIONAL :: kfind                       ! 0  target cell already found
693      !                                                                    ! 1  target to be determined
694      REAL(wp),     INTENT(in   )           :: pdvol, pdsal, pdtem, pratio ! vol/sal/tem increment
695      !                                                                    ! and ratio in case increment span over multiple cells.
696      !!----------------------------------------------------------------------
697      INTEGER :: ifind
698      !!----------------------------------------------------------------------
699      !
700      ! increment position
701      kpts = kpts + 1
702      !
703      ! define if we need to look for closest valid wet cell (no neighbours or neigbourg on halo)
704      IF ( PRESENT(kfind) ) THEN
705         ifind = kfind
706      ELSE
707         ifind = ( 1 - tmask_h(ki,kj) ) * tmask(ki,kj,kk)
708      END IF
709      !
710      ! update isfpts structure
711      sisfpts(kpts) = isfcons(mig(ki), mjg(kj), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, glamt(ki,kj), gphit(ki,kj), ifind )
712      !
713   END SUBROUTINE update_isfpts
714   !
715   SUBROUTINE get_correction( ki, kj, kk, plon, plat, pvolinc, psalinc, pteminc, kfind)
716      !!---------------------------------------------------------------------
717      !!                  ***  ROUTINE get_correction  ***
718      !!
719      !! ** Action : - Find the closest valid cell if needed (wet and not on the halo)
720      !!             - Scale the correction depending of pratio (case where multiple wet neigbourgs)
721      !!             - Fill the correction array
722      !!
723      !!----------------------------------------------------------------------
724      INTEGER , INTENT(in) :: ki, kj, kk, kfind        ! target point indices
725      REAL(wp), INTENT(in) :: plon, plat               ! target point lon/lat
726      REAL(wp), INTENT(in) :: pvolinc, pteminc,psalinc ! correction increment for vol/temp/salt
727      !!----------------------------------------------------------------------
728      INTEGER :: jj, ji, iig, ijg
729      !!----------------------------------------------------------------------
730      !
731      ! define global indice of correction location
732      iig = ki ; ijg = kj
733      IF ( kfind == 1 ) CALL dom_ngb( plon, plat, iig, ijg,'T', kk)
734      !
735      ! fill the correction array
736      DO jj = mj0(ijg),mj1(ijg)
737         DO ji = mi0(iig),mi1(iig)
738            ! correct the vol_flx and corresponding heat/salt flx in the closest cell
739            risfcpl_cons_vol(ji,jj,kk)        =  risfcpl_cons_vol(ji,jj,kk       ) + pvolinc
740            risfcpl_cons_tsc(ji,jj,kk,jp_sal) =  risfcpl_cons_tsc(ji,jj,kk,jp_sal) + psalinc
741            risfcpl_cons_tsc(ji,jj,kk,jp_tem) =  risfcpl_cons_tsc(ji,jj,kk,jp_tem) + pteminc
742         END DO
743      END DO
744
745   END SUBROUTINE get_correction
746
747END MODULE isfcpl
Note: See TracBrowser for help on using the repository browser.