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 branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90 @ 8135

Last change on this file since 8135 was 8135, checked in by timgraham, 7 years ago

Added changes to code in NST_SRC from dev_r5803_UKMO_AGRIF_Vert_interp
  • 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   PUBLIC reconstructandremap
20
21   !                                              !!* Namelist namagrif: AGRIF parameters
22   LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !:
23   INTEGER , PUBLIC ::   nn_cln_update = 3         !: update frequency
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
29   !                                              !!! OLD namelist names
30   INTEGER , PUBLIC ::   nbcline = 0               !: update counter
31   INTEGER , PUBLIC ::   nbclineupdate             !: update frequency
32   REAL(wp), PUBLIC ::   visc_tra                  !: sponge coeff. for tracers
33   REAL(wp), PUBLIC ::   visc_dyn                  !: sponge coeff. for dynamics
34
35   LOGICAL , PUBLIC :: spongedoneT = .FALSE.       !: tracer   sponge layer indicator
36   LOGICAL , PUBLIC :: spongedoneU = .FALSE.       !: dynamics sponge layer indicator
37   LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE.     !: if true: first step
38   LOGICAL , PUBLIC :: lk_agrif_doupd = .TRUE.     !: if true: send update from current grid
39   LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE.    !: if true: print debugging info
40
41   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_tsn
42# if defined key_top
43   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_trn
44# endif
45   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_u
46   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_v
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,  DIMENSION(:,:) :: fsaht_spu, fsaht_spv !: sponge diffusivities
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,  DIMENSION(:,:) :: fsahm_spt, fsahm_spf !: sponge viscosities
49
50   ! Barotropic arrays used to store open boundary data during
51   ! time-splitting loop:
52   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::  ubdy_w, vbdy_w, hbdy_w
53   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::  ubdy_e, vbdy_e, hbdy_e
54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::  ubdy_n, vbdy_n, hbdy_n
55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::  ubdy_s, vbdy_s, hbdy_s
56
57   INTEGER :: tsn_id                                                  ! AGRIF profile for tracers interpolation and update
58   INTEGER :: un_interp_id, vn_interp_id                              ! AGRIF profiles for interpolations
59   INTEGER :: un_update_id, vn_update_id                              ! AGRIF profiles for udpates
60   INTEGER :: tsn_sponge_id, un_sponge_id, vn_sponge_id               ! AGRIF profiles for sponge layers
61# if defined key_top
62   INTEGER :: trn_id, trn_sponge_id
63# endif 
64   INTEGER :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id
65   INTEGER :: ub2b_update_id, vb2b_update_id
66   INTEGER :: e3t_id, e1u_id, e2v_id, sshn_id
67   INTEGER :: scales_t_id
68# if defined key_zdftke
69   INTEGER :: avt_id, avm_id, en_id
70# endif 
71   INTEGER :: umsk_id, vmsk_id
72   INTEGER :: kindic_agr
73
74   !!----------------------------------------------------------------------
75   !! NEMO/NST 3.3.1 , NEMO Consortium (2011)
76   !! $Id$
77   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
78   !!----------------------------------------------------------------------
79CONTAINS
80
81   INTEGER FUNCTION agrif_oce_alloc()
82      !!----------------------------------------------------------------------
83      !!                ***  FUNCTION agrif_oce_alloc  ***
84      !!----------------------------------------------------------------------
85      INTEGER, DIMENSION(2) :: ierr
86      !!----------------------------------------------------------------------
87      ierr(:) = 0
88      !
89      ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj),   &
90         &      fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj),   &
91         &      tabspongedone_tsn(jpi,jpj),           &
92# if defined key_top         
93         &      tabspongedone_trn(jpi,jpj),           &
94# endif         
95         &      tabspongedone_u  (jpi,jpj),           &
96         &      tabspongedone_v  (jpi,jpj), STAT = ierr(1) )
97
98      ALLOCATE( ubdy_w(jpj), vbdy_w(jpj), hbdy_w(jpj),   &
99         &      ubdy_e(jpj), vbdy_e(jpj), hbdy_e(jpj),   & 
100         &      ubdy_n(jpi), vbdy_n(jpi), hbdy_n(jpi),   & 
101         &      ubdy_s(jpi), vbdy_s(jpi), hbdy_s(jpi), STAT = ierr(2) )
102
103      agrif_oce_alloc = MAXVAL(ierr)
104      !
105   END FUNCTION agrif_oce_alloc
106
107      subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)     
108      implicit none
109      integer N, Nout
110      real tabin(N), tabout(Nout)
111      real hin(N), hout(Nout)
112      real coeffremap(N,3),zwork(N,3)
113      real zwork2(N+1,3)
114      integer k
115      double precision, parameter :: dsmll=1.0d-8 
116      real q,q01,q02,q001,q002,q0
117      real z_win(1:N+1), z_wout(1:Nout+1)
118      real,parameter :: dpthin = 1.D-3
119      integer :: k1, kbox, ktop, ka, kbot
120      real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
121
122      z_win(1)=0.; z_wout(1)= 0.
123      do k=1,N
124       z_win(k+1)=z_win(k)+hin(k)
125      enddo 
126     
127      do k=1,Nout
128       z_wout(k+1)=z_wout(k)+hout(k)       
129      enddo       
130
131        do k=2,N
132          zwork(k,1)=1./(hin(k-1)+hin(k))
133        enddo
134       
135        do k=2,N-1
136          q0 = 1./(hin(k-1)+hin(k)+hin(k+1))
137          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1)
138          zwork(k,3)=q0
139        enddo       
140     
141        do k= 2,N
142        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1))
143        enddo
144       
145        coeffremap(:,1) = tabin(:)
146 
147         do k=2,N-1
148        q001 = hin(k)*zwork2(k+1,1)
149        q002 = hin(k)*zwork2(k,1)       
150        if (q001*q002 < 0) then
151          q001 = 0.
152          q002 = 0.
153        endif
154        q=zwork(k,2)
155        q01=q*zwork2(k+1,1)
156        q02=q*zwork2(k,1)
157        if (abs(q001) > abs(q02)) q001 = q02
158        if (abs(q002) > abs(q01)) q002 = q01
159       
160        q=(q001-q002)*zwork(k,3)
161        q001=q001-q*hin(k+1)
162        q002=q002+q*hin(k-1)
163       
164        coeffremap(k,3)=coeffremap(k,1)+q001
165        coeffremap(k,2)=coeffremap(k,1)-q002
166       
167        zwork2(k,1)=(2.*q001-q002)**2
168        zwork2(k,2)=(2.*q002-q001)**2
169        enddo
170       
171        do k=1,N
172        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then
173        coeffremap(k,3) = coeffremap(k,1)
174        coeffremap(k,2) = coeffremap(k,1)
175        zwork2(k,1) = 0.
176        zwork2(k,2) = 0.
177        endif
178        enddo
179       
180        do k=2,N
181        q002=max(zwork2(k-1,2),dsmll)
182        q001=max(zwork2(k,1),dsmll)
183        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,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 k=1,N
190        q01=zwork2(k+1,3)-coeffremap(k,1)
191        q02=coeffremap(k,1)-zwork2(k,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(k,2)=coeffremap(k,1)-q02
203        coeffremap(k,3)=coeffremap(k,1)+q01
204        enddo
205
206      zbot=0.0
207      kbot=1
208      do k=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(k+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= k+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
244
245            if     (z_wout(k)  .ne.z_win(kbox)   .or.z_wout(k+1).ne.z_win(kbox+1)     ) then
246
247              if     (hin(kbox).gt.dpthin) then
248                q001 = (zbox-z_win(kbox))/hin(kbox)
249                q002 = (ztop-z_win(kbox))/hin(kbox)
250                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
251                q02=q01-1.+(q001+q002)
252                q0=1.-q01-q02
253              else
254                q0 = 1.0
255                q01 = 0.
256                q02 = 0.
257              endif
258          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
259             
260            else
261            tabout(k) = tabin(kbox)
262             
263            endif
264
265          else
266
267            if     (ktop.le.k .and. kbox.ge.k) then
268              ka = k
269            elseif (kbox-ktop.ge.3) then
270              ka = (kbox+ktop)/2
271            elseif (hin(ktop).ge.hin(kbox)) then
272              ka = ktop
273            else
274              ka = kbox
275            endif !choose ka
276
277            offset=coeffremap(ka,1)
278
279            qtop = z_win(ktop+1)-ztop !partial layer thickness
280            if     (hin(ktop).gt.dpthin) then
281              q=(ztop-z_win(ktop))/hin(ktop)
282              q01=q*(q-1.)
283              q02=q01+q
284              q0=1-q01-q02           
285            else
286              q0 = 1.
287              q01 = 0.
288              q02 = 0.
289            endif
290           
291            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
292
293            do k1= ktop+1,kbox-1
294              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
295            enddo !k1
296
297           
298            qbot = zbox-z_win(kbox) !partial layer thickness
299            if     (hin(kbox).gt.dpthin) then
300              q=qbot/hin(kbox)
301              q01=(q-1.)**2
302              q02=q01-1.+q
303              q0=1-q01-q02                           
304            else
305              q0 = 1.0
306              q01 = 0.
307              q02 = 0.
308            endif
309           
310            tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
311           
312            rpsum=1.0d0/zthk
313              tabout(k)=offset+tsum*rpsum
314             
315          endif !single or multiple layers
316        else
317        if (k==1) then
318        write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1)
319        endif
320         tabout(k) = tabout(k-1)
321         
322        endif !normal:thin layer
323      enddo !k
324           
325      return
326      end subroutine reconstructandremap
327     
328#endif
329   !!======================================================================
330END MODULE agrif_oce
Note: See TracBrowser for help on using the repository browser.