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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90 @ 11574

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

#2222, import changes from dev_r10973_AGRIF-01_jchanut_small_jpi_jpj (i.e. #2199)

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