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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce.F90 @ 11047

Last change on this file since 11047 was 11047, checked in by acc, 5 years ago

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