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

Last change on this file since 10023 was 10023, checked in by gm, 5 years ago

#1911 (ENHANCE-04): RK3 branch - step II.2 bug correction in dynnxt + domvvl_RK3 creation

File size: 17.8 KB
Line 
1MODULE iscplrst
2   !!======================================================================
3   !!                       ***  MODULE  iscplrst  ***
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   !!----------------------------------------------------------------------
10   !!   iscpl_stp          : step management
11   !!   iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
12   !!----------------------------------------------------------------------
13   USE oce             ! global tra/dyn variable
14   USE dom_oce         ! ocean space and time domain
15   USE domvvl          !
16   USE domwri          ! ocean space and time domain
17   USE phycst          ! physical constants
18   USE sbc_oce         ! surface boundary condition variables
19   USE iscplini        ! ice sheet coupling: initialisation
20   USE iscplhsb        ! ice sheet coupling: conservation
21   !
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
26   USE lbclnk          ! communication
27
28   IMPLICIT NONE
29   PRIVATE
30   
31   PUBLIC   iscpl_stp          ! step management
32   !!
33   !! * Substitutions 
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
37   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
38   !! Software governed by the CeCILL licence     (./LICENSE)
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, ssh(Nnn) (wetting/drying)   
48      !!              compute correction term if needed
49      !!
50      !!----------------------------------------------------------------------
51      INTEGER  ::   inum0
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
56      CHARACTER(20) :: cfile
57      !!----------------------------------------------------------------------
58      !
59      !                       ! get restart variable
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)
68      !
69      CALL iscpl_init()       ! read namelist
70      !                       ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
71      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
72      !
73      IF ( ln_hsb ) THEN      ! compute correction if conservation needed
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         
78      !                       ! create  a domain file
79      IF( ln_meshmask .AND. ln_iscpl )   CALL dom_wri
80      !
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
90      !
91      l_1st_euler = .TRUE.    ! next step is an euler time step
92      !
93      !                       ! set _b and _n variables equal
94      tsb(:,:,:,:) = tsn(:,:,:,:)
95      ub (:,:,:)   = un (:,:,:)
96      vb (:,:,:)   = vn (:,:,:)
97      ssh(:,:,Nbb) = ssh(:,:,Nnn)
98      !
99      !                       ! set _b and _n vertical scale factor equal
100      e3t_b (:,:,:) = e3t_n (:,:,:)
101      e3u_b (:,:,:) = e3u_n (:,:,:)
102      e3v_b (:,:,:) = e3v_n (:,:,:)
103      !
104      e3uw_b (:,:,:) = e3uw_n (:,:,:)
105      e3vw_b (:,:,:) = e3vw_n (:,:,:)
106      gdept_b(:,:,:) = gdept_n(:,:,:)
107      gdepw_b(:,:,:) = gdepw_n(:,:,:)
108      hu_b   (:,:)   = hu_n   (:,:)
109      hv_b   (:,:)   = hv_n   (:,:)
110      r1_hu_b(:,:)   = r1_hu_n(:,:)
111      r1_hv_b(:,:)   = r1_hv_n(:,:)
112      !
113   END SUBROUTINE iscpl_stp
114
115
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 ssh(Nnn) 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
128      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
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      !!
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
140      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
141      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
142      REAL(wp), DIMENSION(jpi,jpj) ::   z_ssh_h0, zsshu, zsshv, zsshf
143      !!----------------------------------------------------------------------
144      !
145      !                 ! mask value to be sure
146      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
147      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
148      !
149      !                 ! compute wmask
150      zwmaskn(:,:,1) = tmask   (:,:,1)
151      zwmaskb(:,:,1) = ptmask_b(:,:,1)
152      DO jk = 2, jpk
153         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
154         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
155      END DO
156      !   
157      !                 ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
158      ssh  (:,:,Nbb) = ssh(:,:,Nnn)
159      zssh0(:,:)     = ssh(:,:,Nnn)
160      zsmask0(:,:)   = psmask_b(:,:)
161      zsmask1(:,:)   = psmask_b(:,:)
162      DO iz = 1, 10                 ! need to be tuned (configuration dependent) (OK for ISOMIP+)
163         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
164         DO jj = 2,jpj-1
165            DO ji = fs_2, fs_jpim1   ! vector opt.
166               summsk = zsmask0(ji+1,jj)+zsmask0(ji-1,jj)+zsmask0(ji,jj+1)+zsmask0(ji,jj-1)
167               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
168                  ssh(ji,jj,Nnn)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj)     &
169                  &              + zssh0(ji-1,jj)*zsmask0(ji-1,jj)     &
170                  &              + zssh0(ji,jj+1)*zsmask0(ji,jj+1)     &
171                  &              + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk
172                  zsmask1(ji,jj) = 1._wp
173               ENDIF
174            END DO
175         END DO
176         CALL lbc_lnk_multi( ssh(:,:,Nnn), 'T', 1., zsmask1, 'T', 1. )
177         zssh0(:,:) = ssh(:,:,Nnn)
178         zsmask0 = zsmask1
179      END DO
180      ssh(:,:,Nnn) = ssh(:,:,Nnn) * ssmask(:,:)
181
182!!gm BUGs....
183!
184!     for me ht_0, hu_0, hv_0 and hf_0 should be recomputed whatever the value of ln_linssh
185!     further more mask at all grid-point should be recomputed
186!     and mikt, u, v, f also...
187!
188!     perhaps, not, if dom_cfg.nc file has been modified....
189!
190!     Pierre we should discuss of that !
191
192
193
194
195!=============================================================================
196!PM: Is this needed since introduction of VVL by default?
197      IF ( .NOT.ln_linssh ) THEN
198         ! Reconstruction of all vertical scale factors at now time steps
199         ! ======================================================================
200
201
202         CALL ctl_stop( 'iscplrst : gm: here there is a BUG: not all required fields are defined')
203
204         !
205         !                          !* NOW fields :
206         CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
207         !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
208         !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
209         !                                !      gdept, gdepw, gde3w         (from 1 to jpk  )
210         !
211
212
213     
214!!gm Question : bug ????
215      ! at this stage  is ht_0, r1_ht0 already computed ????
216      ! I have the feeling that this should be done in a different manner...
217      ! that is iscpl_rst_interpol  should be call directly in restart !!!
218     
219      ! in my changes here I assume that ht_0 AND  r1_ht_0  have been updated
220      ! Note that the former calculation were using ht_0  so if it as not been updated ===>>> BUG
221!!gm
222      ENDIF
223
224!=============================================================================
225! compute velocity
226! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
227      ub(:,:,:) = un(:,:,:)
228      vb(:,:,:) = vn(:,:,:)
229      DO jk = 1, jpk
230         DO jj = 1, jpj
231            DO ji = 1, jpi
232               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);
233               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);
234            END DO
235         END DO
236      END DO
237
238! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
239! compute barotropic velocity now and after
240      zbub(:,:)   = SUM( ub(:,:,:)*pe3u_b(:,:,:) , DIM=3 )
241      zbvb(:,:)   = SUM( vb(:,:,:)*pe3v_b(:,:,:) , DIM=3 )
242      zbun(:,:)   = SUM( un(:,:,:)* e3u_n(:,:,:) , DIM=3 )
243      zbvn(:,:)   = SUM( vn(:,:,:)* e3v_n(:,:,:) , DIM=3 )
244
245      ! new water column
246      zhu1=0.0_wp ;
247      zhv1=0.0_wp ;
248      DO jk  = 1,jpk
249        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
250        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
251      END DO
252     
253      ! compute correction     
254      zucorr = 0._wp
255      zvcorr = 0._wp
256      DO jj = 1, jpj
257         DO ji = 1, jpi
258            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
259               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
260            END IF
261            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
262               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
263            END IF
264         END DO
265      END DO 
266     
267      ! update velocity
268      DO jk  = 1, jpk
269         un(:,:,jk) = ( un(:,:,jk) - zucorr(:,:) ) * umask(:,:,jk)
270         vn(:,:,jk) = ( vn(:,:,jk) - zvcorr(:,:) ) * vmask(:,:,jk)
271      END DO
272
273!=============================================================================
274      ! compute temp and salt
275      ! compute new tn and sn if we open a new cell
276      tsb  (:,:,:,:) = tsn   (:,:,:,:)
277      zts0 (:,:,:,:) = tsn   (:,:,:,:)
278      ztmask1(:,:,:) = ptmask_b(:,:,:)
279      ztmask0(:,:,:) = ptmask_b(:,:,:)
280      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
281         DO jk = 1,jpk-1
282            zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
283            DO jj = 2,jpj-1
284               DO ji = fs_2,fs_jpim1
285                  jip1=ji+1; jim1=ji-1;
286                  jjp1=jj+1; jjm1=jj-1;
287                  summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
288                  IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
289                     !! horizontal basic extrapolation
290                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
291                         &            +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
292                         &            +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
293                         &            +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
294                     tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
295                         &            +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
296                         &            +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
297                         &            +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
298                     ztmask1(ji,jj,jk)=1
299                  ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
300                     !! vertical extrapolation if horizontal extrapolation failed
301                     jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
302                     summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
303                     IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
304                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
305                            &            +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / summsk
306                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
307                            &            +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / summsk
308                        ztmask1(ji,jj,jk)=1._wp
309                     ENDIF
310                  ENDIF
311               END DO
312            END DO
313         END DO
314         !
315         CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
316         !
317         ! update
318         zts0(:,:,:,:) = tsn(:,:,:,:)
319         ztmask0 = ztmask1
320         !
321      END DO
322
323      ! mask new tsn field
324      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
325      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
326
327      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
328      !PM: Is this IF needed since change to VVL by default
329      IF (.NOT.ln_linssh) THEN
330         DO jk = 2,jpk-1
331            DO jj = 1,jpj
332               DO ji = 1,jpi
333                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND.   &
334                  &      (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
335                     !compute weight
336                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
337                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
338                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  ))
339                     IF (zdz .LT. 0._wp) THEN
340                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
341                     END IF
342                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
343                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
344                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk)
345                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
346                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
347                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk)
348                  END IF
349               END DO
350            END DO
351         END DO               
352      END IF
353
354      ! closed pool
355      ! -----------------------------------------------------------------------------------------
356      ! case we open a cell but no neigbour cells available to get an estimate of T and S
357      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
358         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
359         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
360         umask(:,:,:) = 0._wp
361         vmask(:,:,:) = 0._wp
362      END WHERE
363     
364      ! set mbkt and mikt to 1 in thiese location
365      WHERE (SUM(tmask,dim=3) == 0)
366         mbkt(:,:) = 1   ;   mbku(:,:)=1   ;   mbkv(:,:) = 1
367         mikt(:,:) = 1   ;   miku(:,:)=1   ;   mikv(:,:) = 1
368      END WHERE
369      ! -------------------------------------------------------------------------------------------
370      ! compute new tn and sn if we close cell
371      ! nothing to do
372      !
373   END SUBROUTINE iscpl_rst_interpol
374
375   !!======================================================================
376END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.