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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplhsb.F90 @ 10978

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

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : DOM and sshwzv.F90. (sshwzv.F90 will be rewritten somewhat when we rewrite the time-level swapping but I've done it anyway). Passes all non-AGRIF SETTE tests.
  • Property svn:keywords set to Id
File size: 16.1 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( Kbb, Kmm, 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*e3tn-tsb*e3tb and e3tn-e3tb
51      !!----------------------------------------------------------------------
52      INTEGER                     , INTENT(in ) :: Kbb, Kmm    !! time level indices
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      REAL(wp) ::   summsk, zsum , zsumn, zjip1_ratio  , zjim1_ratio, zdtem, zde3t, z1_rdtiscpl
63      REAL(wp) ::   zarea , zsum1, zsumb, zjjp1_ratio  , zjjm1_ratio, zdsal
64      REAL(wp), DIMENSION(jpi,jpj)        ::   zdssh   ! workspace
65      REAL(wp), DIMENSION(:), ALLOCATABLE ::   zlon, zlat
66      REAL(wp), DIMENSION(:), ALLOCATABLE ::   zcorr_vol, zcorr_tem, zcorr_sal
67      INTEGER , DIMENSION(:), ALLOCATABLE ::   ixpts, iypts, izpts, inpts
68      INTEGER :: jpts, npts
69      !!----------------------------------------------------------------------
70
71      ! get imbalance (volume heat and salt)
72      ! initialisation difference
73      zde3t = 0._wp   ;   zdsal = 0._wp   ;   zdtem = 0._wp
74
75      ! initialisation correction term
76      pvol_flx(:,:,:  ) = 0._wp
77      pts_flx (:,:,:,:) = 0._wp
78     
79      z1_rdtiscpl = 1._wp / prdt_iscpl 
80
81      ! mask ts(:,:,:,:,Kmm) and ts(:,:,:,:,Kbb)
82      ts(:,:,:,jp_tem,Kbb) = ts(:,:,:,jp_tem,Kbb) * ptmask_b(:,:,:)
83      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) *  tmask  (:,:,:)
84      ts(:,:,:,jp_sal,Kbb) = ts(:,:,:,jp_sal,Kbb) * ptmask_b(:,:,:)
85      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) *  tmask  (:,:,:)
86
87      !==============================================================================
88      ! diagnose the heat, salt and volume input and compute the correction variable
89      !==============================================================================
90
91      !
92      zdssh(:,:) = ssh(:,:,Kmm) * ssmask(:,:) - ssh(:,:,Kbb) * psmask_b(:,:)
93      IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
94     
95      DO jk = 1,jpk-1
96         DO jj = 2,jpj-1
97            DO ji = fs_2,fs_jpim1
98               IF (tmask_h(ji,jj) == 1._wp) THEN
99
100                  ! volume differences
101                  zde3t = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
102
103                  ! heat diff
104                  zdtem = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) *  tmask  (ji,jj,jk)   &
105                        - ts(ji,jj,jk,jp_tem,Kbb) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
106                  ! salt diff
107                  zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) *  tmask  (ji,jj,jk)   &
108                        - ts(ji,jj,jk,jp_sal,Kbb) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
109               
110                  ! shh changes
111                  IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN
112                     zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
113                     zdssh(ji,jj) = 0._wp
114                  END IF
115
116                  ! volume, heat and salt differences in each cell
117                  pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t * z1_rdtiscpl
118                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal * z1_rdtiscpl 
119                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem * z1_rdtiscpl
120
121                  ! case where we close a cell: check if the neighbour cells are wet
122                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
123
124                     jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
125
126                     zsum =   e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
127                       &    + e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk)
128
129                     IF ( zsum /= 0._wp ) THEN
130                        zjip1_ratio   = e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
131                        zjim1_ratio   = e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
132                        zjjp1_ratio   = e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
133                        zjjm1_ratio   = e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
134
135                        pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
136                        pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
137                        pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
138                        pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
139                        pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
140                        pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
141                        pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
142                        pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
143                        pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
144                        pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
145                        pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
146                        pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
147
148                        ! set to 0 the cell we distributed over neigbourg cells
149                        pvol_flx(ji,jj,jk       ) = 0._wp
150                        pts_flx (ji,jj,jk,jp_sal) = 0._wp
151                        pts_flx (ji,jj,jk,jp_tem) = 0._wp
152
153                     ELSE IF (zsum == 0._wp ) THEN
154                        ! case where we close a cell and no adjacent cell open
155                        ! check if the cell beneath is wet
156                        IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
157                           pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
158                           pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
159                           pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
160
161                           ! set to 0 the cell we distributed over neigbourg cells
162                           pvol_flx(ji,jj,jk       ) = 0._wp
163                           pts_flx (ji,jj,jk,jp_sal) = 0._wp
164                           pts_flx (ji,jj,jk,jp_tem) = 0._wp
165                        ELSE
166                        ! case no adjacent cell on the horizontal and on the vertical
167                           IF ( lwp ) THEN   ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
168                                             ! open for narea== 1 (lwp=T)
169                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
170                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk
171                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea
172                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
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( Kmm, 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      INTEGER                   , INTENT(in   ) ::   Kmm      ! time level index
311      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
312      !!
313      INTEGER  ::   ji, jj, jk   ! dummy loop indices
314      !!----------------------------------------------------------------------
315      !
316      DO jk = 1, jpk
317         DO jj = 1, jpj
318            DO ji = 1, jpi
319               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t(ji,jj,jk,Kmm)
320            END DO
321         END DO
322      END DO
323      !
324   END SUBROUTINE iscpl_div
325
326END MODULE iscplhsb
Note: See TracBrowser for help on using the repository browser.