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

Last change on this file since 11931 was 11931, checked in by mathiot, 10 months ago

ENHANCE-02_ISF_nemo: add comments, improve memory usage of ln_isfcpl_cons option, fix issue in ISOMIP+ configuration

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