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.
iscplrst.F90 in branches/UKMO/dev_r8600_nn_etau_options/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_r8600_nn_etau_options/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 @ 8875

Last change on this file since 8875 was 8875, checked in by davestorkey, 6 years ago

UKMO/dev_r8600_nn_etau_options branch: remove SVN keywords.

File size: 20.8 KB
RevLine 
[5790]1MODULE iscplrst
2   !!======================================================================
[7646]3   !!                       ***  MODULE  iscplrst  ***
[5790]4   !! Ocean forcing: update the restart file in case of ice sheet/ocean coupling
5   !!=====================================================================
6   !! History :  NEMO  ! 2015-01 P. Mathiot: original
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
[5835]10   !!   iscpl_stp          : step management
[5790]11   !!   iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
12   !!----------------------------------------------------------------------
13   USE dom_oce         ! ocean space and time domain
14   USE domwri          ! ocean space and time domain
15   USE domvvl, ONLY : dom_vvl_interpol
16   USE phycst          ! physical constants
17   USE sbc_oce         ! surface boundary condition variables
18   USE oce             ! global tra/dyn variable
19   USE in_out_manager  ! I/O manager
20   USE iom             ! I/O module
21   USE lib_mpp         ! MPP library
22   USE lib_fortran     ! MPP library
23   USE wrk_nemo        ! Memory allocation
[5823]24   USE lbclnk          ! communication
25   USE iscplini        ! ice sheet coupling: initialisation
26   USE iscplhsb        ! ice sheet coupling: conservation
[5790]27
28   IMPLICIT NONE
29   PRIVATE
30   
[5823]31   PUBLIC   iscpl_stp          ! step management
[5790]32   !!
33   !! * Substitutions 
[5920]34#  include "vectopt_loop_substitute.h90"
[5790]35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[8875]37   !! $Id$
[5790]38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE iscpl_stp
43      !!----------------------------------------------------------------------
44      !!                   ***  ROUTINE iscpl_stp  ***
45      !!
46      !! ** Purpose : compute initialisation
47      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
48      !!              compute correction term if needed
49      !!
50      !!----------------------------------------------------------------------
[5945]51      INTEGER  ::   inum0
[7646]52      REAL(wp), DIMENSION(:,:  ), POINTER ::   zsmask_b
53      REAL(wp), DIMENSION(:,:,:), POINTER ::   ztmask_b, zumask_b, zvmask_b
54      REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3t_b  , ze3u_b  , ze3v_b 
55      REAL(wp), DIMENSION(:,:,:), POINTER ::   zdepw_b
[5945]56      CHARACTER(20) :: cfile
[5823]57      !!----------------------------------------------------------------------
[5790]58
[7646]59      CALL wrk_alloc(jpi,jpj,jpk,   ztmask_b, zumask_b, zvmask_b) ! mask before
60      CALL wrk_alloc(jpi,jpj,jpk,   ze3t_b  , ze3u_b  , ze3v_b  ) ! e3   before
61      CALL wrk_alloc(jpi,jpj,jpk,   zdepw_b )
62      CALL wrk_alloc(jpi,jpj,       zsmask_b                    )
[5790]63
64
65      !! get restart variable
66      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b   ) ! need to extrapolate T/S
67      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b   ) ! need to correct barotropic velocity
68      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b   ) ! need to correct barotropic velocity
69      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b   ) ! need to correct barotropic velocity
[6074]70      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:) )  ! need to compute temperature correction
71      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b(:,:,:) )  ! need to correct barotropic velocity
72      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b(:,:,:) )  ! need to correct barotropic velocity
73      CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:) ) ! need to interpol vertical profile (vvl)
[5790]74
75      !! read namelist
76      CALL iscpl_init()
77
78      !!  ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
79      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
80
81      !! compute correction if conservation needed
82      IF ( ln_hsb ) THEN
83         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
84         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
85      END IF
86         
87      !! print mesh/mask
[7646]88      IF( nn_msh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
[5790]89
90      IF ( ln_hsb ) THEN
91         cfile='correction'
92         cfile = TRIM( cfile )
93         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
94         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
95         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
96         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
97         CALL iom_close ( inum0 )
98      END IF
99
[7646]100      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask_b,zumask_b,zvmask_b ) 
101      CALL wrk_dealloc(jpi,jpj,jpk,   ze3t_b  ,ze3u_b  ,ze3v_b   ) 
102      CALL wrk_dealloc(jpi,jpj,jpk,   zdepw_b                    )
103      CALL wrk_dealloc(jpi,jpj,       zsmask_b                   )
[5790]104
[5802]105      !! next step is an euler time step
106      neuler = 0
[5835]107
[5802]108      !! set _b and _n variables equal
109      tsb (:,:,:,:) = tsn (:,:,:,:)
[7646]110      ub  (:,:,:)   = un  (:,:,:)
111      vb  (:,:,:)   = vn  (:,:,:)
112      sshb(:,:)     = sshn(:,:)
[5835]113
[5802]114      !! set _b and _n vertical scale factor equal
[6074]115      e3t_b (:,:,:) = e3t_n (:,:,:)
116      e3u_b (:,:,:) = e3u_n (:,:,:)
117      e3v_b (:,:,:) = e3v_n (:,:,:)
[5945]118
[7646]119      e3uw_b (:,:,:) = e3uw_n (:,:,:)
120      e3vw_b (:,:,:) = e3vw_n (:,:,:)
121      gdept_b(:,:,:) = gdept_n(:,:,:)
[6080]122      gdepw_b(:,:,:) = gdepw_n(:,:,:)
[7646]123      hu_b   (:,:)   = hu_n   (:,:)
124      hv_b   (:,:)   = hv_n   (:,:)
125      r1_hu_b(:,:)   = r1_hu_n(:,:)
126      r1_hv_b(:,:)   = r1_hv_n(:,:)
[5790]127      !
128   END SUBROUTINE iscpl_stp
[7646]129
130
[5790]131   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
132      !!----------------------------------------------------------------------
133      !!                   ***  ROUTINE iscpl_rst_interpol  ***
134      !!
135      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
136      !!                compute 2d fields of heat, salt and volume correction
137      !!
138      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
139      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
140      !!----------------------------------------------------------------------
141      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
142      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
[5823]143      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
[5790]144      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
145      !!
146      INTEGER :: ji, jj, jk, iz          !! loop index
147      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
148      !!
149      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
150      REAL(wp):: zdz, zdzm1, zdzp1
151      !!
152      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
153      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
[5945]154      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, zhu1, zhv1
[5790]155      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
156      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
157      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
158      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
[5945]159      !!----------------------------------------------------------------------
[5823]160
161      !! allocate variables
[5790]162      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
[5835]163      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        ) 
[5790]164      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
165      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
166      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
167      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
[5945]168      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, zhu1, zhv1             ) 
[5823]169
[5790]170      !! mask value to be sure
171      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
172      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
173     
174      ! compute wmask
175      zwmaskn(:,:,1) = tmask   (:,:,1)
176      zwmaskb(:,:,1) = ptmask_b(:,:,1)
177      DO jk = 2,jpk
178         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
179         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
180      END DO
181           
182      ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
183      sshb (:,:)=sshn(:,:)
184      zssh0(:,:)=sshn(:,:)
[5802]185      zsmask0(:,:) = psmask_b(:,:)
186      zsmask1(:,:) = psmask_b(:,:)
[5823]187      DO iz = 1,10    ! need to be tuned (configuration dependent) (OK for ISOMIP+)
[5790]188         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
[5920]189         DO jj = 2,jpj-1
190            DO ji = fs_2, fs_jpim1   ! vector opt.
[5790]191               jip1=ji+1; jim1=ji-1;
192               jjp1=jj+1; jjm1=jj-1;
193               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
[5945]194               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
[5820]195                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
196                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
197                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
198                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
[5823]199                  zsmask1(ji,jj)=1._wp
[5790]200               END IF
201            END DO
202         END DO
[5823]203         CALL lbc_lnk(sshn,'T',1._wp)
204         CALL lbc_lnk(zsmask1,'T',1._wp)
[5790]205         zssh0   = sshn
206         zsmask0 = zsmask1
207      END DO
208      sshn(:,:) = sshn(:,:) * ssmask(:,:)
209
210!=============================================================================
[6075]211!PM: Is this needed since introduction of VVL by default?
212      IF (.NOT.ln_linssh) THEN
[5820]213      ! Reconstruction of all vertical scale factors at now time steps
214      ! =============================================================================
215      ! Horizontal scale factor interpolations
216      ! --------------------------------------
[5790]217         DO jk = 1,jpk
[5820]218            DO jj=1,jpj
219               DO ji=1,jpi
220                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
[8329]221                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) /   &
222                     &   ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
[5820]223                  ENDIF
224               END DO
225            END DO
[5790]226         END DO
227
[6074]228         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
229         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
230         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
[5790]231
232      ! Vertical scale factor interpolations
233      ! ------------------------------------
[6074]234         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
235         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
236         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
[5790]237
238      ! t- and w- points depth
239      ! ----------------------
[6074]240         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
241         gdepw_n(:,:,1) = 0.0_wp
242         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
[5790]243         DO jj = 1,jpj
244            DO ji = 1,jpi
245               DO jk = 2,mikt(ji,jj)-1
[6074]246                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
247                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
248                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
[5790]249               END DO
[5945]250               IF (mikt(ji,jj) > 1) THEN
[5790]251                  jk = mikt(ji,jj)
[6074]252                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
253                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
254                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
[5790]255               END IF
256               DO jk = mikt(ji,jj)+1, jpk
[6074]257                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
258                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
259                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
[5790]260               END DO
261            END DO
262         END DO
263
[5823]264      ! t-, u- and v- water column thickness
265      ! ------------------------------------
[6075]266         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
[5790]267         DO jk = 1, jpkm1
[6074]268            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
269            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
270            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
[5790]271         END DO
272         !                                        ! Inverse of the local depth
[6074]273         r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
274         r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
[5790]275
276      END IF
277
278!=============================================================================
279! compute velocity
280! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
281      ub(:,:,:)=un(:,:,:)
282      vb(:,:,:)=vn(:,:,:)
283      DO jk = 1,jpk
[5920]284         DO jj = 1,jpj
285            DO ji = 1,jpi
[6074]286               un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk);
287               vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk);
[5790]288            END DO
289         END DO
290      END DO
[5802]291
[5790]292! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
293! compute barotropic velocity now and after
294      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
295      zbub(:,:)   = SUM(ztrp,DIM=3)
296      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
297      zbvb(:,:)   = SUM(ztrp,DIM=3)
[6074]298      ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 
[5790]299      zbun(:,:)   = SUM(ztrp,DIM=3)
[6074]300      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
[5790]301      zbvn(:,:)   = SUM(ztrp,DIM=3)
[5802]302
[5820]303      ! new water column
[5945]304      zhu1=0.0_wp ;
305      zhv1=0.0_wp ;
[5790]306      DO jk  = 1,jpk
[6074]307        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
308        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
[5790]309      END DO
[5820]310     
311      ! compute correction     
[5790]312      zucorr = 0._wp
313      zvcorr = 0._wp
[5920]314      DO jj = 1,jpj
315         DO ji = 1,jpi
[5945]316            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
317               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
[5790]318            END IF
[5945]319            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
320               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
[5790]321            END IF
322         END DO
323      END DO 
[5820]324     
325      ! update velocity
[5790]326      DO jk  = 1,jpk
327         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
328         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
329      END DO
330
331!=============================================================================
332      ! compute temp and salt
333      ! compute new tn and sn if we open a new cell
334      tsb (:,:,:,:) = tsn(:,:,:,:)
335      zts0(:,:,:,:) = tsn(:,:,:,:)
336      ztmask1(:,:,:) = ptmask_b(:,:,:)
337      ztmask0(:,:,:) = ptmask_b(:,:,:)
[5920]338      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
[5790]339          DO jk = 1,jpk-1
340              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
[5920]341              DO jj = 2,jpj-1
342                 DO ji = fs_2,fs_jpim1
[5790]343                      jip1=ji+1; jim1=ji-1;
344                      jjp1=jj+1; jjm1=jj-1;
[5820]345                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
[5945]346                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
[5820]347                      !! horizontal basic extrapolation
348                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
349                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
350                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
351                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
352                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
353                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
354                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
355                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
[5790]356                         ztmask1(ji,jj,jk)=1
[5945]357                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
[5920]358                      !! vertical extrapolation if horizontal extrapolation failed
[5790]359                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
360                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
[5945]361                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
[5790]362                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
363                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
364                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
365                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
[5823]366                            ztmask1(ji,jj,jk)=1._wp
[5790]367                         END IF
368                      END IF
369                  END DO
370              END DO
371          END DO
[5823]372         
373          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
374          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
[5945]375          CALL lbc_lnk(ztmask1,     'T',1._wp)
[5823]376
377          ! update
[5790]378          zts0(:,:,:,:) = tsn(:,:,:,:)
379          ztmask0 = ztmask1
[5823]380
[5790]381      END DO
[5823]382
383      ! mask new tsn field
[5790]384      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
385      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
386
[5802]387      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
[6075]388      !PM: Is this IF needed since change to VVL by default
389      IF (.NOT.ln_linssh) THEN
[5790]390         DO jk = 2,jpk-1
391            DO jj = 1,jpj
392               DO ji = 1,jpi
[8329]393                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND.   &
394                  &      (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
[5790]395                     !compute weight
[6074]396                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
397                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
398                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  ))
[5823]399                     IF (zdz .LT. 0._wp) THEN
[5820]400                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
[5790]401                     END IF
402                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
403                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
[6074]404                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk)
[5790]405                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
406                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
[6074]407                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk)
[5790]408                  END IF
409               END DO
410            END DO
411         END DO               
412      END IF
413
[5823]414      ! closed pool
415      ! -----------------------------------------------------------------------------------------
416      ! case we open a cell but no neigbour cells available to get an estimate of T and S
417      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
[5945]418         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
[5823]419         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
420         umask(:,:,:) = 0._wp
421         vmask(:,:,:) = 0._wp
[5790]422      END WHERE
[5823]423     
424      ! set mbkt and mikt to 1 in thiese location
[5790]425      WHERE (SUM(tmask,dim=3) == 0)
426         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
427         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
428      END WHERE
[5823]429      ! -------------------------------------------------------------------------------------------
[5790]430      ! compute new tn and sn if we close cell
431      ! nothing to do
[5802]432      !
433      ! deallocation tmp arrays
[5790]434      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
435      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
436      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d              ) 
437      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
438      CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
439      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
[5945]440      CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , zhu1 , zhv1          ) 
[7646]441      !
[5790]442   END SUBROUTINE iscpl_rst_interpol
443
[7646]444   !!======================================================================
[5790]445END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.