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 branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/DOM – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/DOM/iscplhsb.F90 @ 8568

Last change on this file since 8568 was 8568, checked in by gm, 7 years ago

#1911 (ENHANCE-09): PART I.2 - _NONE option + remove zts + see associated wiki page
File size: 15.5 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/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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      !
91      zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
92      IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
93     
94      DO jk = 1,jpk-1
95         DO jj = 2,jpj-1
96            DO ji = fs_2,fs_jpim1
97               IF (tmask_h(ji,jj) == 1._wp) THEN
98
99                  ! volume differences
100                  zde3t = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
101
102                  ! heat diff
103                  zdtem = tsn(ji,jj,jk,jp_tem) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
104                        - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
105                  ! salt diff
106                  zdsal = tsn(ji,jj,jk,jp_sal) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
107                        - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
108               
109                  ! shh changes
110                  IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN
111                     zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
112                     zdssh(ji,jj) = 0._wp
113                  END IF
114
115                  ! volume, heat and salt differences in each cell
116                  pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t * z1_rdtiscpl
117                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal * z1_rdtiscpl 
118                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem * z1_rdtiscpl
119
120                  ! case where we close a cell: check if the neighbour cells are wet
121                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
122
123                     jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
124
125                     zsum =   e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
126                       &    + e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk)
127
128                     IF ( zsum /= 0._wp ) THEN
129                        zjip1_ratio   = e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
130                        zjim1_ratio   = e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
131                        zjjp1_ratio   = e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
132                        zjjm1_ratio   = e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
133
134                        pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
135                        pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
136                        pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
137                        pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
138                        pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
139                        pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
140                        pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
141                        pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
142                        pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
143                        pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
144                        pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
145                        pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
146
147                        ! set to 0 the cell we distributed over neigbourg cells
148                        pvol_flx(ji,jj,jk       ) = 0._wp
149                        pts_flx (ji,jj,jk,jp_sal) = 0._wp
150                        pts_flx (ji,jj,jk,jp_tem) = 0._wp
151
152                     ELSE IF (zsum == 0._wp ) THEN
153                        ! case where we close a cell and no adjacent cell open
154                        ! check if the cell beneath is wet
155                        IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
156                           pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
157                           pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
158                           pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
159
160                           ! set to 0 the cell we distributed over neigbourg cells
161                           pvol_flx(ji,jj,jk       ) = 0._wp
162                           pts_flx (ji,jj,jk,jp_sal) = 0._wp
163                           pts_flx (ji,jj,jk,jp_tem) = 0._wp
164                        ELSE
165                        ! case no adjacent cell on the horizontal and on the vertical
166                           IF ( lwp ) THEN   ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
167                                             ! open for narea== 1 (lwp=T)
168                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
169                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk
170                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea
171                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
172                           ENDIF
173                        ! We deal with these points later.
174                        END IF
175                     END IF
176                  END IF
177               END IF
178            END DO
179         END DO
180      END DO
181
182      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
183      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
184      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
185
186      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point
187      ! allocation and initialisation of the list of problematic point
188      ALLOCATE( inpts(jpnij) )
189      inpts(:) = 0
190
191      ! fill narea location with the number of problematic point
192      DO jk = 1,jpk-1
193         DO jj = 2,jpj-1
194            DO ji = fs_2,fs_jpim1
195               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
196                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
197                  inpts(narea) = inpts(narea) + 1 
198               END IF
199            END DO
200         END DO
201      END DO
202
203      ! build array of total problematic point on each cpu (share to each cpu)
204      CALL mpp_max(inpts,jpnij) 
205
206      ! size of the new variable
207      npts  = SUM(inpts)   
208     
209      ! allocation of the coordinates, correction, index vector for the problematic points
210      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
211      ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20_wp ; zlat(:) = -1.0e20_wp
212      zcorr_vol(:) = -1.0e20_wp
213      zcorr_sal(:) = -1.0e20_wp
214      zcorr_tem(:) = -1.0e20_wp
215
216      ! fill new variable
217      jpts = SUM(inpts(1:narea-1))
218      DO jk = 1,jpk-1
219         DO jj = 2,jpj-1
220            DO ji = fs_2,fs_jpim1
221               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
222                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
223                  jpts = jpts + 1  ! positioning in the inpts vector for the area narea
224                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
225                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
226                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
227                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
228                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
229
230                  ! set flx to 0 (safer)
231                  pvol_flx(ji,jj,jk       ) = 0.0_wp
232                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
233                  pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
234               END IF
235            END DO
236         END DO
237      END DO
238
239      ! build array of total problematic point on each cpu (share to each cpu)
240      ! point coordinates
241      CALL mpp_max(zlat ,npts)
242      CALL mpp_max(zlon ,npts)
243      CALL mpp_max(izpts,npts)
244
245      ! correction values
246      CALL mpp_max(zcorr_vol,npts)
247      CALL mpp_max(zcorr_sal,npts)
248      CALL mpp_max(zcorr_tem,npts)
249
250      ! put correction term in the closest cell         
251      DO jpts = 1,npts
252         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
253         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
254            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
255               jk = izpts(jpts)
256
257               IF (tmask_h(ji,jj) == 1._wp) THEN
258                  ! correct the vol_flx in the closest cell
259                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
260                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
261                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
262
263                  ! set correction to 0
264                  zcorr_vol(jpts) = 0.0_wp
265                  zcorr_sal(jpts) = 0.0_wp
266                  zcorr_tem(jpts) = 0.0_wp
267               END IF
268            END DO
269         END DO
270      END DO
271
272      ! deallocate variables
273      DEALLOCATE(inpts)
274      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
275   
276      ! add contribution store on the hallo (lbclnk remove one of the contribution)
277      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
278      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
279      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
280
281      ! compute sum over the halo and set it to 0.
282      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
283      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
284      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
285      !
286   END SUBROUTINE iscpl_cons
287
288
289   SUBROUTINE iscpl_div( phdivn )
290      !!----------------------------------------------------------------------
291      !!                  ***  ROUTINE iscpl_div  ***
292      !!
293      !! ** Purpose :   update the horizontal divergenc
294      !!
295      !! ** Method  :
296      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
297      !!
298      !! ** Action  :   phdivn   increase by the iscpl correction term
299      !!----------------------------------------------------------------------
300      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
301      !!
302      INTEGER  ::   ji, jj, jk   ! dummy loop indices
303      !!----------------------------------------------------------------------
304      !
305      DO jk = 1, jpk
306         DO jj = 1, jpj
307            DO ji = 1, jpi
308               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t_n(ji,jj,jk)
309            END DO
310         END DO
311      END DO
312      !
313   END SUBROUTINE iscpl_div
314
315END MODULE iscplhsb
Note: See TracBrowser for help on using the repository browser.