source: NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST/agrif_oce.F90 @ 11219

Last change on this file since 11219 was 11219, checked in by jchanut, 16 months ago

#2199
1) Define aditionnal arrays to correct the time interpolation of barotropic arrays in corners. Since multiple stages in the time interpolation are necessary, overlapping segments in corners give wrong results otherwise (corrects stage 2 in previous commit)..
2) Added subroutine to correct time extrapolated fluxes at bdy in time splitting routine (updates stage 3 in previous commit).
3) Completly remove non-specified open boundary case. Boundares are now exactly set from parent (no more filtering nor extrapolation for outgoing flows).
At this stage, use of nbondi, nbondj variables has been completly removed.

  • Property svn:keywords set to Id
File size: 11.7 KB
RevLine 
[782]1MODULE agrif_oce
[1605]2   !!======================================================================
[782]3   !!                       ***  MODULE agrif_oce  ***
[1605]4   !! AGRIF :   define in memory AGRIF variables
[782]5   !!----------------------------------------------------------------------
[1605]6   !! History :  2.0  ! 2007-12  (R. Benshila)  Original code
[782]7   !!----------------------------------------------------------------------
[1605]8#if defined key_agrif
[782]9   !!----------------------------------------------------------------------
[1605]10   !!   'key_agrif'                                              AGRIF zoom
11   !!----------------------------------------------------------------------
[782]12   USE par_oce      ! ocean parameters
13   USE dom_oce      ! domain parameters
[5656]14
[782]15   IMPLICIT NONE
[2715]16   PRIVATE
[782]17
[2715]18   PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90
[9031]19#if defined key_vertical
20   PUBLIC reconstructandremap ! remapping routine
21#endif   
[1605]22   !                                              !!* Namelist namagrif: AGRIF parameters
[5656]23   LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !:
[10425]24   INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 2  !: Sponge width (in number of parent grid points)
[5656]25   REAL(wp), PUBLIC ::   rn_sponge_tra = 2800.     !: sponge coeff. for tracers
26   REAL(wp), PUBLIC ::   rn_sponge_dyn = 2800.     !: sponge coeff. for dynamics
27   LOGICAL , PUBLIC ::   ln_chk_bathy  = .FALSE.   !: check of parent bathymetry
[1605]28   !                                              !!! OLD namelist names
29   REAL(wp), PUBLIC ::   visc_tra                  !: sponge coeff. for tracers
30   REAL(wp), PUBLIC ::   visc_dyn                  !: sponge coeff. for dynamics
[782]31
[5656]32   LOGICAL , PUBLIC :: spongedoneT = .FALSE.       !: tracer   sponge layer indicator
33   LOGICAL , PUBLIC :: spongedoneU = .FALSE.       !: dynamics sponge layer indicator
34   LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE.     !: if true: first step
35   LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE.    !: if true: print debugging info
[1605]36
[5656]37   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_tsn
38# if defined key_top
39   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_trn
40# endif
41   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_u
42   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_v
[11219]43   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utint_stage
44   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtint_stage
[9019]45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsaht_spu, fsaht_spv !: sponge diffusivities
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsahm_spt, fsahm_spf !: sponge viscosities
[2715]47
[9019]48   ! Barotropic arrays used to store open boundary data during time-splitting loop:
[11205]49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ubdy, vbdy, hbdy
[5656]50
[9019]51
52   INTEGER, PUBLIC :: tsn_id                                                  ! AGRIF profile for tracers interpolation and update
53   INTEGER, PUBLIC :: un_interp_id, vn_interp_id                              ! AGRIF profiles for interpolations
54   INTEGER, PUBLIC :: un_update_id, vn_update_id                              ! AGRIF profiles for udpates
55   INTEGER, PUBLIC :: tsn_sponge_id, un_sponge_id, vn_sponge_id               ! AGRIF profiles for sponge layers
[5656]56# if defined key_top
[9019]57   INTEGER, PUBLIC :: trn_id, trn_sponge_id
[5656]58# endif 
[9019]59   INTEGER, PUBLIC :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id
60   INTEGER, PUBLIC :: ub2b_update_id, vb2b_update_id
61   INTEGER, PUBLIC :: e3t_id, e1u_id, e2v_id, sshn_id
62   INTEGER, PUBLIC :: scales_t_id
63   INTEGER, PUBLIC :: avt_id, avm_id, en_id                ! TKE related identificators
64   INTEGER, PUBLIC :: umsk_id, vmsk_id
65   INTEGER, PUBLIC :: kindic_agr
66   
[1605]67   !!----------------------------------------------------------------------
[9598]68   !! NEMO/NST 4.0 , NEMO Consortium (2018)
[1605]69   !! $Id$
[10068]70   !! Software governed by the CeCILL license (see ./LICENSE)
[2715]71   !!----------------------------------------------------------------------
72CONTAINS
73
74   INTEGER FUNCTION agrif_oce_alloc()
75      !!----------------------------------------------------------------------
76      !!                ***  FUNCTION agrif_oce_alloc  ***
77      !!----------------------------------------------------------------------
[5656]78      INTEGER, DIMENSION(2) :: ierr
79      !!----------------------------------------------------------------------
80      ierr(:) = 0
81      !
[11219]82      ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj),     &
83         &      fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj),     &
84         &      tabspongedone_tsn(jpi,jpj),                 &
85         &      utint_stage(jpi,jpj), vtint_stage(jpi,jpj), &
[5656]86# if defined key_top         
87         &      tabspongedone_trn(jpi,jpj),           &
88# endif         
89         &      tabspongedone_u  (jpi,jpj),           &
90         &      tabspongedone_v  (jpi,jpj), STAT = ierr(1) )
91
[11205]92      ALLOCATE( ubdy(jpi,jpj), vbdy(jpi,jpj), hbdy(jpi,jpj), STAT = ierr(2) )
[5656]93
94      agrif_oce_alloc = MAXVAL(ierr)
95      !
[2715]96   END FUNCTION agrif_oce_alloc
97
[9031]98#if defined key_vertical
99   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)     
100      !!----------------------------------------------------------------------
101      !!                ***  FUNCTION reconstructandremap  ***
102      !!----------------------------------------------------------------------
103      IMPLICIT NONE
104      INTEGER N, Nout
105      REAL(wp) tabin(N), tabout(Nout)
106      REAL(wp) hin(N), hout(Nout)
107      REAL(wp) coeffremap(N,3),zwork(N,3)
108      REAL(wp) zwork2(N+1,3)
109      INTEGER jk
110      DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8 
111      REAL(wp) q,q01,q02,q001,q002,q0
112      REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)
113      REAL(wp),PARAMETER :: dpthin = 1.D-3
114      INTEGER :: k1, kbox, ktop, ka, kbot
115      REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
116
117      z_win(1)=0.; z_wout(1)= 0.
118      DO jk=1,N
119         z_win(jk+1)=z_win(jk)+hin(jk)
120      ENDDO 
121     
122      DO jk=1,Nout
123         z_wout(jk+1)=z_wout(jk)+hout(jk)       
124      ENDDO       
125
126      DO jk=2,N
127         zwork(jk,1)=1./(hin(jk-1)+hin(jk))
128      ENDDO
129       
130      DO jk=2,N-1
131         q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))
132         zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)
133         zwork(jk,3)=q0
134      ENDDO       
135     
136      DO jk= 2,N
137         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
138      ENDDO
139       
140      coeffremap(:,1) = tabin(:)
141 
142      DO jk=2,N-1
143         q001 = hin(jk)*zwork2(jk+1,1)
144         q002 = hin(jk)*zwork2(jk,1)       
145         IF (q001*q002 < 0) then
146            q001 = 0.
147            q002 = 0.
148         ENDIF
149         q=zwork(jk,2)
150         q01=q*zwork2(jk+1,1)
151         q02=q*zwork2(jk,1)
152         IF (abs(q001) > abs(q02)) q001 = q02
153         IF (abs(q002) > abs(q01)) q002 = q01
154       
155         q=(q001-q002)*zwork(jk,3)
156         q001=q001-q*hin(jk+1)
157         q002=q002+q*hin(jk-1)
158       
159         coeffremap(jk,3)=coeffremap(jk,1)+q001
160         coeffremap(jk,2)=coeffremap(jk,1)-q002
161       
162         zwork2(jk,1)=(2.*q001-q002)**2
163         zwork2(jk,2)=(2.*q002-q001)**2
164      ENDDO
165       
166      DO jk=1,N
167         IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN
168            coeffremap(jk,3) = coeffremap(jk,1)
169            coeffremap(jk,2) = coeffremap(jk,1)
170            zwork2(jk,1) = 0.
171            zwork2(jk,2) = 0.
172         ENDIF
173      ENDDO
174       
175      DO jk=2,N
176         q002=max(zwork2(jk-1,2),dsmll)
177         q001=max(zwork2(jk,1),dsmll)
178         zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)
179      ENDDO
180       
181      zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
182      zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
183 
184      DO jk=1,N
185         q01=zwork2(jk+1,3)-coeffremap(jk,1)
186         q02=coeffremap(jk,1)-zwork2(jk,3)
187         q001=2.*q01
188         q002=2.*q02
189         IF (q01*q02<0) then
190            q01=0.
191            q02=0.
192         ELSEIF (abs(q01)>abs(q002)) then
193            q01=q002
194         ELSEIF (abs(q02)>abs(q001)) then
195            q02=q001
196         ENDIF
197         coeffremap(jk,2)=coeffremap(jk,1)-q02
198         coeffremap(jk,3)=coeffremap(jk,1)+q01
199      ENDDO
200
201      zbot=0.0
202      kbot=1
203      DO jk=1,Nout
204         ztop=zbot  !top is bottom of previous layer
205         ktop=kbot
206         IF     (ztop.GE.z_win(ktop+1)) then
207            ktop=ktop+1
208         ENDIF
209       
210         zbot=z_wout(jk+1)
211         zthk=zbot-ztop
212
213         IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN
214
215            kbot=ktop
216            DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
217               kbot=kbot+1
218            ENDDO
219            zbox=zbot
220            DO k1= jk+1,Nout
221               IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN
222                  exit !thick layer
223               ELSE
224                  zbox=z_wout(k1+1)  !include thin adjacent layers
225                  IF(zbox.EQ.z_wout(Nout+1)) THEN
226                     exit !at bottom
227                  ENDIF
228               ENDIF
229            ENDDO
230            zthk=zbox-ztop
231
232            kbox=ktop
233            DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
234               kbox=kbox+1
235            ENDDO
236         
237            IF(ktop.EQ.kbox) THEN
238               IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN
239                  IF(hin(kbox).GT.dpthin) THEN
240                     q001 = (zbox-z_win(kbox))/hin(kbox)
241                     q002 = (ztop-z_win(kbox))/hin(kbox)
242                     q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
243                     q02=q01-1.+(q001+q002)
244                     q0=1.-q01-q02
245                  ELSE
246                     q0 = 1.0
247                     q01 = 0.
248                     q02 = 0.
249                  ENDIF
250                  tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
251               ELSE
252                  tabout(jk) = tabin(kbox)
253               ENDIF
254            ELSE
255               IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
256                  ka = jk
257               ELSEIF (kbox-ktop.GE.3) THEN
258                  ka = (kbox+ktop)/2
259               ELSEIF (hin(ktop).GE.hin(kbox)) THEN
260                  ka = ktop
261               ELSE
262                  ka = kbox
263               ENDIF !choose ka
264   
265               offset=coeffremap(ka,1)
266   
267               qtop = z_win(ktop+1)-ztop !partial layer thickness
268               IF(hin(ktop).GT.dpthin) THEN
269                  q=(ztop-z_win(ktop))/hin(ktop)
270                  q01=q*(q-1.)
271                  q02=q01+q
272                  q0=1-q01-q02           
273               ELSE
274                  q0 = 1.
275                  q01 = 0.
276                  q02 = 0.
277               ENDIF
278               
279               tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
280   
281               DO k1= ktop+1,kbox-1
282                  tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
283               ENDDO !k1
284               
285               qbot = zbox-z_win(kbox) !partial layer thickness
286               IF(hin(kbox).GT.dpthin) THEN
287                  q=qbot/hin(kbox)
288                  q01=(q-1.)**2
289                  q02=q01-1.+q
290                  q0=1-q01-q02                           
291               ELSE
292                  q0 = 1.0
293                  q01 = 0.
294                  q02 = 0.
295               ENDIF
296             
297               tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
298               
299               rpsum=1.0d0/zthk
300               tabout(jk)=offset+tsum*rpsum
301                 
302            ENDIF !single or multiple layers
303         ELSE
304            IF (jk==1) THEN
305               write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)
306            ENDIF
307            tabout(jk) = tabout(jk-1)
308             
309         ENDIF !normal:thin layer
310      ENDDO !jk
311           
312      return
313      end subroutine reconstructandremap
[2715]314#endif
[9031]315
316#endif
[1605]317   !!======================================================================
[782]318END MODULE agrif_oce
Note: See TracBrowser for help on using the repository browser.