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

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90 @ 9939

Last change on this file since 9939 was 9939, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche

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