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

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

#1911 (ENHANCE-04): RK3 branch - step II.3 remove e3uw_$ e3vw_$, except e3.w_0 and use only after e3 in dyn/trazdf

File size: 17.6 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      gdept_b(:,:,:) = gdept_n(:,:,:)
105      gdepw_b(:,:,:) = gdepw_n(:,:,:)
106      hu_b   (:,:)   = hu_n   (:,:)
107      hv_b   (:,:)   = hv_n   (:,:)
108      r1_hu_b(:,:)   = r1_hu_n(:,:)
109      r1_hv_b(:,:)   = r1_hv_n(:,:)
110      !
111   END SUBROUTINE iscpl_stp
112
113
114   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
115      !!----------------------------------------------------------------------
116      !!                   ***  ROUTINE iscpl_rst_interpol  ***
117      !!
118      !! ** Purpose :   compute new tn, sn, un, vn and ssh(Nnn) in case of evolving geometry of ice shelves
119      !!                compute 2d fields of heat, salt and volume correction
120      !!
121      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
122      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
123      !!----------------------------------------------------------------------
124      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
125      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
126      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
127      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
128      !!
129      INTEGER :: ji, jj, jk, iz          !! loop index
130      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
131      !!
132      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
133      REAL(wp):: zdz, zdzm1, zdzp1
134      !!
135      REAL(wp), DIMENSION(jpi,jpj)          :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t
136      REAL(wp), DIMENSION(jpi,jpj)          :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1
137      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn
138      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
139      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
140      REAL(wp), DIMENSION(jpi,jpj) ::   z_ssh_h0, zsshu, zsshv, zsshf
141      !!----------------------------------------------------------------------
142      !
143      !                 ! mask value to be sure
144      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
145      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
146      !
147      !                 ! compute wmask
148      zwmaskn(:,:,1) = tmask   (:,:,1)
149      zwmaskb(:,:,1) = ptmask_b(:,:,1)
150      DO jk = 2, jpk
151         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
152         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
153      END DO
154      !   
155      !                 ! compute new ssh if we open a full water column (average of the closest neigbourgs) 
156      ssh  (:,:,Nbb) = ssh(:,:,Nnn)
157      zssh0(:,:)     = ssh(:,:,Nnn)
158      zsmask0(:,:)   = psmask_b(:,:)
159      zsmask1(:,:)   = psmask_b(:,:)
160      DO iz = 1, 10                 ! need to be tuned (configuration dependent) (OK for ISOMIP+)
161         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
162         DO jj = 2,jpj-1
163            DO ji = fs_2, fs_jpim1   ! vector opt.
164               summsk = zsmask0(ji+1,jj)+zsmask0(ji-1,jj)+zsmask0(ji,jj+1)+zsmask0(ji,jj-1)
165               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
166                  ssh(ji,jj,Nnn)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj)     &
167                  &              + zssh0(ji-1,jj)*zsmask0(ji-1,jj)     &
168                  &              + zssh0(ji,jj+1)*zsmask0(ji,jj+1)     &
169                  &              + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk
170                  zsmask1(ji,jj) = 1._wp
171               ENDIF
172            END DO
173         END DO
174         CALL lbc_lnk_multi( ssh(:,:,Nnn), 'T', 1., zsmask1, 'T', 1. )
175         zssh0(:,:) = ssh(:,:,Nnn)
176         zsmask0 = zsmask1
177      END DO
178      ssh(:,:,Nnn) = ssh(:,:,Nnn) * ssmask(:,:)
179
180!!gm BUGs....
181!
182!     for me ht_0, hu_0, hv_0 and hf_0 should be recomputed whatever the value of ln_linssh
183!     further more mask at all grid-point should be recomputed
184!     and mikt, u, v, f also...
185!
186!     perhaps, not, if dom_cfg.nc file has been modified....
187!
188!     Pierre we should discuss of that !
189
190
191
192
193!=============================================================================
194!PM: Is this needed since introduction of VVL by default?
195      IF ( .NOT.ln_linssh ) THEN
196         ! Reconstruction of all vertical scale factors at now time steps
197         ! ======================================================================
198
199
200         CALL ctl_stop( 'iscplrst : gm: here there is a BUG: not all required fields are defined')
201
202         !
203         !                          !* NOW fields :
204         CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
205         !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
206         !                                !      e3w, gdept, gdepw, gde3w    (from 1 to jpk  )
207         !
208
209
210     
211!!gm Question : bug ????
212      ! at this stage  is ht_0, r1_ht0 already computed ????
213      ! I have the feeling that this should be done in a different manner...
214      ! that is iscpl_rst_interpol  should be call directly in restart !!!
215     
216      ! in my changes here I assume that ht_0 AND  r1_ht_0  have been updated
217      ! Note that the former calculation were using ht_0  so if it as not been updated ===>>> BUG
218!!gm
219      ENDIF
220
221!=============================================================================
222! compute velocity
223! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
224      ub(:,:,:) = un(:,:,:)
225      vb(:,:,:) = vn(:,:,:)
226      DO jk = 1, jpk
227         DO jj = 1, jpj
228            DO ji = 1, jpi
229               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);
230               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);
231            END DO
232         END DO
233      END DO
234
235! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
236! compute barotropic velocity now and after
237      zbub(:,:)   = SUM( ub(:,:,:)*pe3u_b(:,:,:) , DIM=3 )
238      zbvb(:,:)   = SUM( vb(:,:,:)*pe3v_b(:,:,:) , DIM=3 )
239      zbun(:,:)   = SUM( un(:,:,:)* e3u_n(:,:,:) , DIM=3 )
240      zbvn(:,:)   = SUM( vn(:,:,:)* e3v_n(:,:,:) , DIM=3 )
241
242      ! new water column
243      zhu1=0.0_wp ;
244      zhv1=0.0_wp ;
245      DO jk  = 1,jpk
246        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
247        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
248      END DO
249     
250      ! compute correction     
251      zucorr = 0._wp
252      zvcorr = 0._wp
253      DO jj = 1, jpj
254         DO ji = 1, jpi
255            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
256               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
257            END IF
258            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
259               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
260            END IF
261         END DO
262      END DO 
263     
264      ! update velocity
265      DO jk  = 1, jpk
266         un(:,:,jk) = ( un(:,:,jk) - zucorr(:,:) ) * umask(:,:,jk)
267         vn(:,:,jk) = ( vn(:,:,jk) - zvcorr(:,:) ) * vmask(:,:,jk)
268      END DO
269
270!=============================================================================
271      ! compute temp and salt
272      ! compute new tn and sn if we open a new cell
273      tsb  (:,:,:,:) = tsn   (:,:,:,:)
274      zts0 (:,:,:,:) = tsn   (:,:,:,:)
275      ztmask1(:,:,:) = ptmask_b(:,:,:)
276      ztmask0(:,:,:) = ptmask_b(:,:,:)
277      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
278         DO jk = 1,jpk-1
279            zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
280            DO jj = 2,jpj-1
281               DO ji = fs_2,fs_jpim1
282                  jip1=ji+1; jim1=ji-1;
283                  jjp1=jj+1; jjm1=jj-1;
284                  summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
285                  IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
286                     !! horizontal basic extrapolation
287                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
288                         &            +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
289                         &            +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
290                         &            +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
291                     tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
292                         &            +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
293                         &            +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
294                         &            +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
295                     ztmask1(ji,jj,jk)=1
296                  ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
297                     !! vertical extrapolation if horizontal extrapolation failed
298                     jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
299                     summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
300                     IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
301                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
302                            &            +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / summsk
303                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
304                            &            +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / summsk
305                        ztmask1(ji,jj,jk)=1._wp
306                     ENDIF
307                  ENDIF
308               END DO
309            END DO
310         END DO
311         !
312         CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
313         !
314         ! update
315         zts0(:,:,:,:) = tsn(:,:,:,:)
316         ztmask0 = ztmask1
317         !
318      END DO
319
320      ! mask new tsn field
321      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
322      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
323
324      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
325      !PM: Is this IF needed since change to VVL by default
326      IF (.NOT.ln_linssh) THEN
327         DO jk = 2,jpk-1
328            DO jj = 1,jpj
329               DO ji = 1,jpi
330                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND.   &
331                  &      (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
332                     !compute weight
333                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
334                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
335                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  ))
336                     IF (zdz .LT. 0._wp) THEN
337                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
338                     END IF
339                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
340                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
341                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk)
342                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
343                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
344                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk)
345                  END IF
346               END DO
347            END DO
348         END DO               
349      END IF
350
351      ! closed pool
352      ! -----------------------------------------------------------------------------------------
353      ! case we open a cell but no neigbour cells available to get an estimate of T and S
354      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
355         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
356         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
357         umask(:,:,:) = 0._wp
358         vmask(:,:,:) = 0._wp
359      END WHERE
360     
361      ! set mbkt and mikt to 1 in thiese location
362      WHERE (SUM(tmask,dim=3) == 0)
363         mbkt(:,:) = 1   ;   mbku(:,:)=1   ;   mbkv(:,:) = 1
364         mikt(:,:) = 1   ;   miku(:,:)=1   ;   mikv(:,:) = 1
365      END WHERE
366      ! -------------------------------------------------------------------------------------------
367      ! compute new tn and sn if we close cell
368      ! nothing to do
369      !
370   END SUBROUTINE iscpl_rst_interpol
371
372   !!======================================================================
373END MODULE iscplrst
Note: See TracBrowser for help on using the repository browser.