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

Last change on this file since 11823 was 11823, checked in by mathiot, 12 months ago

rm useless USE statement, option compatibility test + minor changes

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