source: branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90 @ 9630

Last change on this file since 9630 was 9630, checked in by antsia, 2 years ago

Code for divergence correction

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