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 @ 5790

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

ice sheet coupling: add/del files for the restart extrapolation (wetting and drying)

File size: 17.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 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
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  )
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       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
86       IF (lwp) PRINT *, 'total volume correction 0 = ',zsum
87       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
88       IF (lwp) PRINT *, 'total heat correction 0 = ',zsum
89       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
90       IF (lwp) PRINT *, 'total salt correction 0 = ',zsum
91
92       ! mask tsn and tsb (should be useless)
93       tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:);
94       tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:);
95
96       ! diagnose non conservation of heat, salt and volume
97       r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt) 
98       zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
99       IF ( lk_vvl ) zssh0 = 0.0_wp
100       DO jk = 1,jpk-1
101          DO ji = 2,jpi-1
102             DO jj = 2,jpj-1
103                ! volume differences
104                zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk);
105               
106                ! shh changes
107                IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN
108                   zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj)
109                   zssh0(ji,jj) = 0._wp
110                END IF
111
112                ! ocean cell now
113                ! case where we open, enlarge or thin a cell :
114                pvol_flx(ji,jj,jk)       =                          zde3t(ji,jj) * r1_tiscpl
115                pts_flx (ji,jj,jk,jp_sal)=   tsn(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl 
116                pts_flx (ji,jj,jk,jp_tem)=   tsn(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl
117             END DO
118          END DO
119       END DO
120       ! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0
121       PRINT *, 'test ', narea, SUM(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt
122       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
123       IF (lwp) PRINT *, 'total volume correction 1 = ',zsum
124       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
125       IF (lwp) PRINT *, 'total heat correction 1 = ',zsum
126       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
127       IF (lwp) PRINT *, 'total salt correction 1 = ',zsum
128
129       zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
130       IF ( lk_vvl ) zssh0 = 0.0_wp
131       DO jk = 1,jpk-1
132          DO ji = 2,jpi-1
133             DO jj = 2,jpj-1
134                ! volume differences
135                zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk);
136               
137                ! shh changes
138                IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN
139                   zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj)
140                   zssh0(ji,jj) = 0._wp
141                END IF
142
143                ! ocean cell before and mask cell now
144                IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
145                   ! case where we close a cell and adjacent cell open
146                   pvol_flx(ji,jj,jk)       =                         zde3t(ji,jj) * r1_tiscpl
147                   pts_flx (ji,jj,jk,jp_sal)=  tsb(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl 
148                   pts_flx (ji,jj,jk,jp_tem)=  tsb(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl
149
150                   jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
151
152                   zsum =   e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
153                     &    + e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) 
154
155                   IF ( zsum .NE. 0._wp ) THEN
156                      zjip1_ratio = e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
157                      zjim1_ratio = e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
158                      zjjp1_ratio = e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
159                      zjjm1_ratio = e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
160
161                      pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
162                      pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
163                      pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
164                      pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
165                      pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
166                      pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
167                      pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
168                      pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
169                      pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
170                      pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
171                      pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
172                      pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
173
174                      ! set to 0 the cell we distributed over neigbourg cells
175                      pvol_flx(ji,jj,jk       ) = 0._wp
176                      pts_flx (ji,jj,jk,jp_sal) = 0._wp
177                      pts_flx (ji,jj,jk,jp_tem) = 0._wp
178
179                   ELSE IF (zsum .EQ. 0._wp ) THEN
180                      ! case where we close a cell and no adjacent cell open
181                      ! check if the cell beneath is wet
182                      IF ( tmask(ji,jj,jk+1) .EQ. 1._wp ) THEN
183                         pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
184                         pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
185                         pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
186
187                         ! set to 0 the cell we distributed over neigbourg cells
188                         pvol_flx(ji,jj,jk       ) = 0._wp
189                         pts_flx (ji,jj,jk,jp_sal) = 0._wp
190                         pts_flx (ji,jj,jk,jp_tem) = 0._wp
191                      ELSE
192                      ! case no adjacent cell on the horizontal and on the vertical
193                         PRINT *,        'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
194                         PRINT *,        '                     ',mig(ji),' ',mjg(jj),' ',jk
195                         PRINT *,        '                     ',ji,' ',jj,' ',jk,' ',narea
196                         PRINT *,        ' we are now looking for the closest wet cell on the horizontal '
197                      ! We deal with this points later.
198                      END IF
199                   END IF
200                END IF
201             END DO
202          END DO
203       END DO
204
205       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
206       IF (lwp) PRINT *, 'total volume correction 2 = ',zsum
207       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
208       IF (lwp) PRINT *, 'total heat correction 2 = ',zsum
209       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
210       IF (lwp) PRINT *, 'total salt correction 2 = ',zsum
211
212       ! allocation and initialisation of the list of problematic point
213       ALLOCATE(vnpts(jpnij))
214       vnpts(:)=0
215
216       ! fill narea location with the number of problematic point
217       DO jk = 1,jpk-1
218          DO ji = 2,jpi-1
219             DO jj = 2,jpj-1
220                IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 &
221                &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 ) THEN
222                   vnpts(narea) = vnpts(narea) + 1 
223                END IF
224             END DO
225          END DO
226       END DO
227
228       ! build array of total problematic point on each cpu (share to each cpu)
229       CALL mpp_max(vnpts,jpnij) 
230
231       ! size of the new variable
232       npts  = SUM(vnpts)   
233       
234       ! allocation of the coordinates, correction, index vector for the problematic points
235       ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
236       ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20
237       zcorr_vol(:) = 0.0_wp
238       zcorr_sal(:) = 0.0_wp
239       zcorr_tem(:) = 0.0_wp
240
241       ! fill new variable
242       jpts = SUM(vnpts(1:narea-1))
243       DO jk = 1,jpk-1
244          DO ji = 2,jpi-1
245             DO jj = 2,jpj-1
246                IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 &
247                &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 ) THEN
248                   jpts = jpts + 1  ! positioning in the vnpts vector for the area narea
249                   PRINT *, 'corrected point ', narea, ji, jj, jk, jpts
250                   ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
251                   zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
252                   zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
253                   zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
254                   zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
255                   ! set flx to 0 (safer)
256                   pvol_flx(ji,jj,jk       ) = 0.0_wp
257                   pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
258                   pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
259                   PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt
260                END IF
261             END DO
262          END DO
263       END DO
264
265       ! build array of total problematic point on each cpu (share to each cpu)
266       CALL mpp_max(zlat ,npts)
267       CALL mpp_max(zlon ,npts)
268       CALL mpp_max(izpts,npts)
269
270       ! put correction term in the closest cell         
271       PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts
272       DO jpts = 1,npts
273          CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
274          PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts)
275          DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
276             DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
277                jk = izpts(jpts)
278                pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
279                pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
280                pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
281             END DO
282          END DO
283       END DO
284       ! deallocate variables
285       DEALLOCATE(vnpts)
286       DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
287     
288       ! add contribution store on the hallo (lbclnk remove one of the contribution)
289       pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
290       pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
291       pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
292
293       CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
294       CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
295       CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
296
297       ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!!
298       zsumn = glob_sum     ( fse3t_n(:,:,:) * tmask  (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt
299       ztmp3d(:,:,:) = 0.0
300       ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
301       zsumb = glob_sum_full(ztmp3d) 
302       zsum  = glob_sum     ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt)
303       IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum
304       ! CHECK salt
305       zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
306       zsumb = glob_sum( tsb(:,:,:,jp_sal) *  pe3t_b(:,:,:) * ptmask_b(:,:,:))
307       zsum  = glob_sum( pts_flx(:,:,:,jp_sal)*rn_fiscpl * rn_rdt)
308       IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum
309       ! CHECK heat
310       zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
311       zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * ptmask_b(:,:,:))
312       zsum  = glob_sum( pts_flx(:,:,:,jp_tem)*rn_fiscpl * rn_rdt)
313       IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum
314    !!
315    CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d ) 
316    CALL wrk_dealloc(jpi,jpj,       zde3t  ) 
317    CALL wrk_dealloc(jpi,jpj,       zssh0  ) 
318   END SUBROUTINE iscpl_cons
319
320   SUBROUTINE iscpl_div( phdivn )
321      !!----------------------------------------------------------------------
322      !!                  ***  ROUTINE iscpl_div  ***
323      !!
324      !! ** Purpose :   update the horizontal divergenc
325      !!
326      !! ** Method  :
327      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
328      !!
329      !! ** Action  :   phdivn   increase by the iscpl correction term
330      !!----------------------------------------------------------------------
331      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
332      !!
333      INTEGER  ::   ji, jj, jk   ! dummy loop indices
334      !!----------------------------------------------------------------------
335      !
336      DO jk = 1, jpk
337         DO jj = 1, jpj
338            DO ji = 1, jpi
339               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk)
340            END DO
341         END DO
342      END DO
343      !
344   END SUBROUTINE iscpl_div
345
346END MODULE iscplhsb
Note: See TracBrowser for help on using the repository browser.