New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
isfcpl.F90 in NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF – NEMO

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

Last change on this file since 11852 was 11852, checked in by mathiot, 4 years ago

ENHANCE-02_ISF_nemo: fix WED025 restartability, finish removing useless USE, remove useless lbc_lnk

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