source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfcpl.F90 @ 12624

Last change on this file since 12624 was 12624, checked in by techene, 8 months ago

all: add e3 substitute and limit precompiled files lines to about 130 character, change key_LF into key_QCO, change module name (dynatfQCO, traatfQCO, stepLF)

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