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/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DOM – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DOM/iscplhsb.F90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago

GMED 462 add flush
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$
38   !! Software governed by the CeCILL license (see ./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      !
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                           IF(lflush) CALL FLUSH(numout)
173                           ENDIF
174                        ! We deal with these points later.
175                        END IF
176                     END IF
177                  END IF
178               END IF
179            END DO
180         END DO
181      END DO
182
183!!gm  ERROR !!!!
184!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
185!
186!      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
187!      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
188!      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
189      STOP ' iscpl_cons:   please modify this module !'
190!!gm end
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('iscplhsb', 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('iscplhsb', zlat ,npts)
247      CALL mpp_max('iscplhsb', zlon ,npts)
248      CALL mpp_max('iscplhsb', izpts,npts)
249
250      ! correction values
251      CALL mpp_max('iscplhsb', zcorr_vol,npts)
252      CALL mpp_max('iscplhsb', zcorr_sal,npts)
253      CALL mpp_max('iscplhsb', 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!!gm  ERROR !!!!
287!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
288!
289!      ! compute sum over the halo and set it to 0.
290!      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
291!      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
292!      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
293!!gm end
294
295      !
296   END SUBROUTINE iscpl_cons
297
298
299   SUBROUTINE iscpl_div( phdivn )
300      !!----------------------------------------------------------------------
301      !!                  ***  ROUTINE iscpl_div  ***
302      !!
303      !! ** Purpose :   update the horizontal divergenc
304      !!
305      !! ** Method  :
306      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
307      !!
308      !! ** Action  :   phdivn   increase by the iscpl correction term
309      !!----------------------------------------------------------------------
310      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
311      !!
312      INTEGER  ::   ji, jj, jk   ! dummy loop indices
313      !!----------------------------------------------------------------------
314      !
315      DO jk = 1, jpk
316         DO jj = 1, jpj
317            DO ji = 1, jpi
318               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t_n(ji,jj,jk)
319            END DO
320         END DO
321      END DO
322      !
323   END SUBROUTINE iscpl_div
324
325END MODULE iscplhsb
Note: See TracBrowser for help on using the repository browser.