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.
iscplhsb.F90 in NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplhsb.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh
File size: 15.9 KB
Line 
1MODULE iscplhsb
2   !!======================================================================
3   !!                       ***  MODULE  iscplhsb  ***
4   !! Ocean forcing: ice sheet/ocean coupling (conservation)
5   !!=====================================================================
6   !! History :  NEMO  ! 2015-01 P. Mathiot: original
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   iscpl_alloc    : variable allocation
11   !!   iscpl_hsb      : compute and store the input of heat/salt/volume
12   !!                    into the system due to the coupling process
13   !!   iscpl_div      : correction of divergence to keep volume conservation
14   !!----------------------------------------------------------------------
15   USE oce             ! global tra/dyn variable
16   USE dom_oce         ! ocean space and time domain
17   USE domwri          ! ocean space and time domain
18   USE domngb          !
19   USE phycst          ! physical constants
20   USE sbc_oce         ! surface boundary condition variables
21   USE iscplini        !
22   !
23   USE in_out_manager  ! I/O manager
24   USE lib_mpp         ! MPP library
25   USE lib_fortran     ! MPP library
26   USE lbclnk          !
27
28   IMPLICIT NONE
29   PRIVATE
30   
31   PUBLIC   iscpl_div   
32   PUBLIC   iscpl_cons       
33   !! * Substitutions 
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
37   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
38   !! Software governed by the CeCILL licence     (./LICENSE)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE iscpl_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl)
43      !!----------------------------------------------------------------------
44      !!                   ***  ROUTINE iscpl_cons  ***
45      !!
46      !! ** Purpose :   compute input into the system during the coupling step
47      !!                compute the correction term
48      !!                compute where the correction have to be applied
49      !!
50      !! ** Method  :   compute tsn*e3t-tsb*e3tb and e3t-e3t_b
51      !!----------------------------------------------------------------------
52      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b    !! mask before
53      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b      !! scale factor before
54      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b    !! mask before
55      REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx     !! corrective flux to have tracer conservation
56      REAL(wp), DIMENSION(:,:,:  ), INTENT(out) :: pvol_flx    !! corrective flux to have volume conservation
57      REAL(wp),                     INTENT(in ) :: prdt_iscpl  !! coupling period
58      !
59      INTEGER  ::   ji  , jj  , jk           ! loop index
60      INTEGER  ::   jip1, jim1, jjp1, jjm1
61      REAL(wp) ::   summsk, zsum , zsumn, zjip1_ratio  , zjim1_ratio, zdtem, zde3t, z1_rdtiscpl
62      REAL(wp) ::   zarea , zsum1, zsumb, zjjp1_ratio  , zjjm1_ratio, zdsal
63      REAL(wp), DIMENSION(jpi,jpj)        ::   zdssh   ! workspace
64      REAL(wp), DIMENSION(:), ALLOCATABLE ::   zlon, zlat
65      REAL(wp), DIMENSION(:), ALLOCATABLE ::   zcorr_vol, zcorr_tem, zcorr_sal
66      INTEGER , DIMENSION(:), ALLOCATABLE ::   ixpts, iypts, izpts, inpts
67      INTEGER :: jpts, npts
68      !!----------------------------------------------------------------------
69
70      ! get imbalance (volume heat and salt)
71      ! initialisation difference
72      zde3t = 0._wp   ;   zdsal = 0._wp   ;   zdtem = 0._wp
73
74      ! initialisation correction term
75      pvol_flx(:,:,:  ) = 0._wp
76      pts_flx (:,:,:,:) = 0._wp
77     
78      z1_rdtiscpl = 1._wp / prdt_iscpl 
79
80      ! mask tsn and tsb
81      tsb(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) * ptmask_b(:,:,:)
82      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) *  tmask  (:,:,:)
83      tsb(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) * ptmask_b(:,:,:)
84      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) *  tmask  (:,:,:)
85
86      !==============================================================================
87      ! diagnose the heat, salt and volume input and compute the correction variable
88      !==============================================================================
89      !
90      zdssh(:,:) = ssh(:,:,Nnn) * ssmask(:,:) - ssh(:,:,Nbb) * psmask_b(:,:)
91      IF (.NOT. ln_linssh )   zdssh = 0._wp     ! already included in the levels by definition
92      !
93      DO jk = 1,jpk-1
94         DO jj = 2,jpj-1
95            DO ji = fs_2,fs_jpim1
96               IF ( tmask_h(ji,jj) == 1._wp ) THEN
97                  !
98                  !                                   ! volume differences
99                  zde3t = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
100
101                  !                                   ! heat diff
102                  zdtem = tsn(ji,jj,jk,jp_tem) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
103                        - tsb(ji,jj,jk,jp_tem) * pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
104                  !                                   ! salt diff
105                  zdsal = tsn(ji,jj,jk,jp_sal) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
106                        - tsb(ji,jj,jk,jp_sal) * pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
107                  !                                   ! shh changes
108                  IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN
109                     zde3t = zde3t + zdssh(ji,jj)     ! zdssh = 0 if vvl
110                     zdssh(ji,jj) = 0._wp
111                  END IF
112
113                  ! volume, heat and salt differences in each cell
114                  pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t * z1_rdtiscpl
115                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal * z1_rdtiscpl 
116                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem * z1_rdtiscpl
117
118                  ! case where we close a cell: check if the neighbour cells are wet
119                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
120
121                     jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
122
123                     zsum =   e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
124                       &    + e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk)
125
126                     IF ( zsum /= 0._wp ) THEN
127                        zjip1_ratio   = e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
128                        zjim1_ratio   = e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
129                        zjjp1_ratio   = e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
130                        zjjm1_ratio   = e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
131
132                        pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
133                        pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
134                        pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
135                        pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
136                        pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
137                        pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
138                        pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
139                        pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
140                        pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
141                        pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
142                        pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
143                        pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
144
145                        ! set to 0 the cell we distributed over neigbourg cells
146                        pvol_flx(ji,jj,jk       ) = 0._wp
147                        pts_flx (ji,jj,jk,jp_sal) = 0._wp
148                        pts_flx (ji,jj,jk,jp_tem) = 0._wp
149
150                     ELSE IF (zsum == 0._wp ) THEN
151                        ! case where we close a cell and no adjacent cell open
152                        ! check if the cell beneath is wet
153                        IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
154                           pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
155                           pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
156                           pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
157
158                           ! set to 0 the cell we distributed over neigbourg cells
159                           pvol_flx(ji,jj,jk       ) = 0._wp
160                           pts_flx (ji,jj,jk,jp_sal) = 0._wp
161                           pts_flx (ji,jj,jk,jp_tem) = 0._wp
162                        ELSE
163                        ! case no adjacent cell on the horizontal and on the vertical
164                           IF ( lwp ) THEN   ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
165                                             ! open for narea== 1 (lwp=T)
166                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
167                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk
168                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea
169                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
170                           ENDIF
171                        ! We deal with these points later.
172                        END IF
173                     END IF
174                  END IF
175               END IF
176            END DO
177         END DO
178      END DO
179
180!!gm  ERROR !!!!
181!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
182!
183!      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
184!      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
185!      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
186      STOP ' iscpl_cons:   please modify this module !'
187!!gm end
188      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point
189      ! allocation and initialisation of the list of problematic point
190      ALLOCATE( inpts(jpnij) )
191      inpts(:) = 0
192
193      ! fill narea location with the number of problematic point
194      DO jk = 1,jpk-1
195         DO jj = 2,jpj-1
196            DO ji = fs_2,fs_jpim1
197               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
198                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
199                  inpts(narea) = inpts(narea) + 1 
200               END IF
201            END DO
202         END DO
203      END DO
204
205      ! build array of total problematic point on each cpu (share to each cpu)
206      CALL mpp_max(inpts,jpnij) 
207
208      ! size of the new variable
209      npts  = SUM(inpts)   
210     
211      ! allocation of the coordinates, correction, index vector for the problematic points
212      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
213      ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20_wp ; zlat(:) = -1.0e20_wp
214      zcorr_vol(:) = -1.0e20_wp
215      zcorr_sal(:) = -1.0e20_wp
216      zcorr_tem(:) = -1.0e20_wp
217
218      ! fill new variable
219      jpts = SUM(inpts(1:narea-1))
220      DO jk = 1,jpk-1
221         DO jj = 2,jpj-1
222            DO ji = fs_2,fs_jpim1
223               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
224                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
225                  jpts = jpts + 1  ! positioning in the inpts vector for the area narea
226                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
227                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
228                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
229                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
230                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
231
232                  ! set flx to 0 (safer)
233                  pvol_flx(ji,jj,jk       ) = 0.0_wp
234                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
235                  pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
236               END IF
237            END DO
238         END DO
239      END DO
240
241      ! build array of total problematic point on each cpu (share to each cpu)
242      ! point coordinates
243      CALL mpp_max(zlat ,npts)
244      CALL mpp_max(zlon ,npts)
245      CALL mpp_max(izpts,npts)
246
247      ! correction values
248      CALL mpp_max(zcorr_vol,npts)
249      CALL mpp_max(zcorr_sal,npts)
250      CALL mpp_max(zcorr_tem,npts)
251
252      ! put correction term in the closest cell         
253      DO jpts = 1,npts
254         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
255         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
256            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
257               jk = izpts(jpts)
258
259               IF (tmask_h(ji,jj) == 1._wp) THEN
260                  ! correct the vol_flx in the closest cell
261                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
262                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
263                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
264
265                  ! set correction to 0
266                  zcorr_vol(jpts) = 0.0_wp
267                  zcorr_sal(jpts) = 0.0_wp
268                  zcorr_tem(jpts) = 0.0_wp
269               END IF
270            END DO
271         END DO
272      END DO
273
274      ! deallocate variables
275      DEALLOCATE(inpts)
276      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
277   
278      ! add contribution store on the hallo (lbclnk remove one of the contribution)
279      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
280      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
281      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
282
283!!gm  ERROR !!!!
284!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
285!
286!      ! compute sum over the halo and set it to 0.
287!      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
288!      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
289!      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
290!!gm end
291
292      !
293   END SUBROUTINE iscpl_cons
294
295
296   SUBROUTINE iscpl_div( phdivn )
297      !!----------------------------------------------------------------------
298      !!                  ***  ROUTINE iscpl_div  ***
299      !!
300      !! ** Purpose :   update the horizontal divergenc
301      !!
302      !! ** Method  :
303      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
304      !!
305      !! ** Action  :   phdivn   increase by the iscpl correction term
306      !!----------------------------------------------------------------------
307      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
308      !!
309      INTEGER  ::   ji, jj, jk   ! dummy loop indices
310      !!----------------------------------------------------------------------
311      !
312      DO jk = 1, jpk
313         DO jj = 1, jpj
314            DO ji = 1, jpi
315               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t_n(ji,jj,jk)
316            END DO
317         END DO
318      END DO
319      !
320   END SUBROUTINE iscpl_div
321
322END MODULE iscplhsb
Note: See TracBrowser for help on using the repository browser.