source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcpl.F90 @ 11553

Last change on this file since 11553 was 11553, checked in by mathiot, 13 months ago

ENHANCE-02_ISF: fix coupling issue (ticket #2142)

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