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/dev_r11943_MERGE_2019/src/OCE/ISF – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isfcpl.F90 @ 12343

Last change on this file since 12343 was 12343, checked in by acc, 5 years ago

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