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_rewrite_time_filterswap/src/NST – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce.F90 @ 11463

Last change on this file since 11463 was 11053, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Merge in latest changes from main branch and finish conversion of "h" variables. NB. This version still doesn't work!
  • Property svn:keywords set to Id
File size: 12.3 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
[9058]28   LOGICAL , PUBLIC ::   lk_agrif_clp  = .FALSE.   !: Force clamped bcs
[1605]29   !                                              !!! OLD namelist names
30   REAL(wp), PUBLIC ::   visc_tra                  !: sponge coeff. for tracers
31   REAL(wp), PUBLIC ::   visc_dyn                  !: sponge coeff. for dynamics
[782]32
[5656]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
[1605]37
[5656]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
[9019]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
[2715]46
[9019]47   ! Barotropic arrays used to store open boundary data during time-splitting loop:
[9116]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
[11053]52   INTEGER , PUBLIC,              SAVE                 ::  Kbb_a, Kmm_a, Krhs_a   !: AGRIF module-specific copies of time-level indices
[5656]53
[9019]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
[5656]59# if defined key_top
[9019]60   INTEGER, PUBLIC :: trn_id, trn_sponge_id
[5656]61# endif 
[9019]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   
[1605]70   !!----------------------------------------------------------------------
[9598]71   !! NEMO/NST 4.0 , NEMO Consortium (2018)
[1605]72   !! $Id$
[10068]73   !! Software governed by the CeCILL license (see ./LICENSE)
[2715]74   !!----------------------------------------------------------------------
75CONTAINS
76
77   INTEGER FUNCTION agrif_oce_alloc()
78      !!----------------------------------------------------------------------
79      !!                ***  FUNCTION agrif_oce_alloc  ***
80      !!----------------------------------------------------------------------
[5656]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
[9116]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) )
[5656]98
99      agrif_oce_alloc = MAXVAL(ierr)
100      !
[2715]101   END FUNCTION agrif_oce_alloc
102
[9031]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
[2715]319#endif
[9031]320
321#endif
[1605]322   !!======================================================================
[782]323END MODULE agrif_oce
Note: See TracBrowser for help on using the repository browser.