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.
agrif_oce.F90 in NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST – NEMO

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

Last change on this file since 11244 was 11244, checked in by jchanut, 5 years ago

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