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

source: NEMO/branches/UKMO/NEMO_4.0_mirror/src/OCE/DOM/iscplhsb.F90 @ 10888

Last change on this file since 10888 was 10888, checked in by davestorkey, 5 years ago

branches/UKMO/NEMO_4.0_mirror : clear SVN keywords

File size: 15.8 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                           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!!gm  ERROR !!!!
183!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
184!
185!      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
186!      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
187!      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
188      STOP ' iscpl_cons:   please modify this module !'
189!!gm end
190      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point
191      ! allocation and initialisation of the list of problematic point
192      ALLOCATE( inpts(jpnij) )
193      inpts(:) = 0
194
195      ! fill narea location with the number of problematic point
196      DO jk = 1,jpk-1
197         DO jj = 2,jpj-1
198            DO ji = fs_2,fs_jpim1
199               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
200                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
201                  inpts(narea) = inpts(narea) + 1 
202               END IF
203            END DO
204         END DO
205      END DO
206
207      ! build array of total problematic point on each cpu (share to each cpu)
208      CALL mpp_max('iscplhsb', inpts,jpnij) 
209
210      ! size of the new variable
211      npts  = SUM(inpts)   
212     
213      ! allocation of the coordinates, correction, index vector for the problematic points
214      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
215      ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20_wp ; zlat(:) = -1.0e20_wp
216      zcorr_vol(:) = -1.0e20_wp
217      zcorr_sal(:) = -1.0e20_wp
218      zcorr_tem(:) = -1.0e20_wp
219
220      ! fill new variable
221      jpts = SUM(inpts(1:narea-1))
222      DO jk = 1,jpk-1
223         DO jj = 2,jpj-1
224            DO ji = fs_2,fs_jpim1
225               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
226                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
227                  jpts = jpts + 1  ! positioning in the inpts vector for the area narea
228                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
229                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
230                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
231                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
232                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
233
234                  ! set flx to 0 (safer)
235                  pvol_flx(ji,jj,jk       ) = 0.0_wp
236                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
237                  pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
238               END IF
239            END DO
240         END DO
241      END DO
242
243      ! build array of total problematic point on each cpu (share to each cpu)
244      ! point coordinates
245      CALL mpp_max('iscplhsb', zlat ,npts)
246      CALL mpp_max('iscplhsb', zlon ,npts)
247      CALL mpp_max('iscplhsb', izpts,npts)
248
249      ! correction values
250      CALL mpp_max('iscplhsb', zcorr_vol,npts)
251      CALL mpp_max('iscplhsb', zcorr_sal,npts)
252      CALL mpp_max('iscplhsb', zcorr_tem,npts)
253
254      ! put correction term in the closest cell         
255      DO jpts = 1,npts
256         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
257         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
258            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
259               jk = izpts(jpts)
260
261               IF (tmask_h(ji,jj) == 1._wp) THEN
262                  ! correct the vol_flx in the closest cell
263                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
264                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
265                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
266
267                  ! set correction to 0
268                  zcorr_vol(jpts) = 0.0_wp
269                  zcorr_sal(jpts) = 0.0_wp
270                  zcorr_tem(jpts) = 0.0_wp
271               END IF
272            END DO
273         END DO
274      END DO
275
276      ! deallocate variables
277      DEALLOCATE(inpts)
278      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
279   
280      ! add contribution store on the hallo (lbclnk remove one of the contribution)
281      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
282      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
283      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
284
285!!gm  ERROR !!!!
286!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
287!
288!      ! compute sum over the halo and set it to 0.
289!      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
290!      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
291!      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
292!!gm end
293
294      !
295   END SUBROUTINE iscpl_cons
296
297
298   SUBROUTINE iscpl_div( phdivn )
299      !!----------------------------------------------------------------------
300      !!                  ***  ROUTINE iscpl_div  ***
301      !!
302      !! ** Purpose :   update the horizontal divergenc
303      !!
304      !! ** Method  :
305      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
306      !!
307      !! ** Action  :   phdivn   increase by the iscpl correction term
308      !!----------------------------------------------------------------------
309      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
310      !!
311      INTEGER  ::   ji, jj, jk   ! dummy loop indices
312      !!----------------------------------------------------------------------
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) / e3t_n(ji,jj,jk)
318            END DO
319         END DO
320      END DO
321      !
322   END SUBROUTINE iscpl_div
323
324END MODULE iscplhsb
Note: See TracBrowser for help on using the repository browser.