source: NEMO/trunk/src/NST/agrif_oce.F90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 2 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks

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