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/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90 @ 5823

Last change on this file since 5823 was 5823, checked in by mathiot, 8 years ago

ice sheet coupling : add comments, rm PRINT statements, cosmetic changes
File size: 15.3 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   !!----------------------------------------------------------------------
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      !!
62      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
63      REAL(wp):: r1_tiscpl
64      REAL(wp):: zjip1_ratio  , zjim1_ratio  , zjjp1_ratio  , zjjm1_ratio
65      !!
66      REAL(wp), DIMENSION(:,:    ), POINTER :: zde3t, zdtem, zdsal
67      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0 
68      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmp3d
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, vnpts
73      INTEGER :: jpts, npts
74
75      CALL wrk_alloc(jpi,jpj,jpk,   ztmp3d )
76      CALL wrk_alloc(jpi,jpj,       zde3t , zdtem, zdsal )
77      CALL wrk_alloc(jpi,jpj,       zssh0  )
78
79    ! get unbalance (volume heat and salt)
80    ! initialisation
81      zde3t   (:,:)     = 0.0_wp
82      pvol_flx(:,:,:  ) = 0.0_wp
83      pts_flx (:,:,:,:) = 0.0_wp
84
85      ! mask tsn and tsb
86      tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:);
87      tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:);
88
89      ! diagnose non conservation of heat, salt and volume
90      r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt) 
91
92      zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
93      IF ( lk_vvl ) zssh0 = 0.0_wp ! already include in the levels by definition
94
95      DO jk = 1,jpk-1
96         DO ji = 2,jpi-1
97            DO jj = 2,jpj-1
98               IF (tmask_h(ji,jj) == 1._wp) THEN
99
100                  ! volume differences
101                  zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
102
103                  ! heat diff
104                  zdtem(ji,jj) = tsn(ji,jj,jk,jp_tem) * fse3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
105                               - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
106                  ! salt diff
107                  zdsal(ji,jj) = tsn(ji,jj,jk,jp_sal) * fse3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
108                               - tsb(ji,jj,jk,jp_sal) * 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(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) ! zssh0 = 0 if vvl
113                     zssh0(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(ji,jj) * r1_tiscpl
118                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal(ji,jj) * r1_tiscpl 
119                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem(ji,jj) * r1_tiscpl
120
121                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
122                     ! case where we close a cell: check if the neighbour cells are wet
123
124                     jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
125
126                     zsum =   e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
127                       &    + e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e12t(jip1,jj  ) * tmask(jip1,jj  ,jk)
128
129                     IF ( zsum .NE. 0._wp ) THEN
130                        zjip1_ratio   = e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
131                        zjim1_ratio   = e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
132                        zjjp1_ratio   = e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
133                        zjjm1_ratio   = e12t(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 .EQ. 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) .EQ. 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                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
168                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk
169                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea
170                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
171                        ! We deal with these points later.
172                        END IF
173                     END IF
174                  END IF
175               END IF
176            END DO
177         END DO
178      END DO
179
180      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
181      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
182      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
183
184      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point
185      ! allocation and initialisation of the list of problematic point
186      ALLOCATE(vnpts(jpnij))
187      vnpts(:)=0
188
189      ! fill narea location with the number of problematic point
190      DO jk = 1,jpk-1
191         DO ji = 2,jpi-1
192            DO jj = 2,jpj-1
193               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
194                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
195                  vnpts(narea) = vnpts(narea) + 1 
196               END IF
197            END DO
198         END DO
199      END DO
200
201      ! build array of total problematic point on each cpu (share to each cpu)
202      CALL mpp_max(vnpts,jpnij) 
203
204      ! size of the new variable
205      npts  = SUM(vnpts)   
206     
207      ! allocation of the coordinates, correction, index vector for the problematic points
208      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
209      ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20
210      zcorr_vol(:) = -1.0e20
211      zcorr_sal(:) = -1.0e20
212      zcorr_tem(:) = -1.0e20
213
214      ! fill new variable
215      jpts = SUM(vnpts(1:narea-1))
216      DO jk = 1,jpk-1
217         DO ji = 2,jpi-1
218            DO jj = 2,jpj-1
219               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
220                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
221                  jpts = jpts + 1  ! positioning in the vnpts vector for the area narea
222                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
223                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
224                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
225                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
226                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
227
228                  ! set flx to 0 (safer)
229                  pvol_flx(ji,jj,jk       ) = 0.0_wp
230                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
231                  pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
232               END IF
233            END DO
234         END DO
235      END DO
236
237      ! build array of total problematic point on each cpu (share to each cpu)
238      ! point coordinates
239      CALL mpp_max(zlat ,npts)
240      CALL mpp_max(zlon ,npts)
241      CALL mpp_max(izpts,npts)
242
243      ! correction values
244      CALL mpp_max(zcorr_vol,npts)
245      CALL mpp_max(zcorr_sal,npts)
246      CALL mpp_max(zcorr_tem,npts)
247
248      ! put correction term in the closest cell         
249      DO jpts = 1,npts
250         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
251         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
252            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
253               jk = izpts(jpts)
254
255               IF (tmask_h(ji,jj) == 1._wp) THEN
256                  ! correct the vol_flx in the closest cell
257                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
258                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
259                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
260
261                  ! set correction to 0
262                  zcorr_vol(jpts) = 0.0_wp
263                  zcorr_sal(jpts) = 0.0_wp
264                  zcorr_tem(jpts) = 0.0_wp
265               END IF
266            END DO
267         END DO
268      END DO
269
270      ! deallocate variables
271      DEALLOCATE(vnpts)
272      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
273   
274      ! add contribution store on the hallo (lbclnk remove one of the contribution)
275      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
276      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
277      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
278
279      ! compute sum over the halo and set it to 0.
280      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
281      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
282      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
283
284      ! deallocate variables
285      CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d ) 
286      CALL wrk_dealloc(jpi,jpj,       zde3t  ) 
287      CALL wrk_dealloc(jpi,jpj,       zssh0  ) 
288   END SUBROUTINE iscpl_cons
289
290   SUBROUTINE iscpl_div( phdivn )
291      !!----------------------------------------------------------------------
292      !!                  ***  ROUTINE iscpl_div  ***
293      !!
294      !! ** Purpose :   update the horizontal divergenc
295      !!
296      !! ** Method  :
297      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
298      !!
299      !! ** Action  :   phdivn   increase by the iscpl correction term
300      !!----------------------------------------------------------------------
301      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
302      !!
303      INTEGER  ::   ji, jj, jk   ! dummy loop indices
304      !!----------------------------------------------------------------------
305      !
306      DO jk = 1, jpk
307         DO jj = 1, jpj
308            DO ji = 1, jpi
309               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk)
310            END DO
311         END DO
312      END DO
313      !
314   END SUBROUTINE iscpl_div
315
316END MODULE iscplhsb
Note: See TracBrowser for help on using the repository browser.