source: trunk/SOURCES/dragging_param_beta_sedim_mod.f90 @ 193

Last change on this file since 193 was 193, checked in by aquiquet, 6 years ago

Correcting a missing variable initialisation: toblim for parameterised beta

File size: 8.5 KB
Line 
1!> \file dragging_neff_contmaj_mod.f90
2!! Defini les zones de stream avec 3 criteres
3!<
4
5!> \namespace  dragging_neff_contmaj
6!! Defini les zones de stream avec 3 criteres
7!! \author ...
8!! \date ...
9!! @note Les trois criteres sont:
10!! @note   * un critere sur la pression effective
11!! @note   * un critere de continuite depuis la cote
12!! @note   * un critere sur la courbure du socle (si negatif, vallees)
13!! @note Used module
14!! @note   - use module3D_phy
15!! @note   - use param_phy_mod
16!<
17
18module dragging_param_beta_sedim
19
20! Defini les zones de stream avec :
21! * un critere sur la pression effective
22! * un critere de continuite depuis la cote
23! * un critere sur la courbure du socle (si negatif, vallees)
24
25use module3d_phy
26use param_phy_mod
27use interface_input
28use io_netcdf_grisli
29use fake_beta_iter_vitbil_mod
30
31implicit none
32
33logical,dimension(nx,ny) :: cote
34
35real,dimension(nx,ny) :: beta_param ! local variable
36
37real :: betamin   ! betamin : (Pa) frottement mini sous les streams
38
39real :: beta_slope        ! = A in:     B = A x Neff ** K
40real :: beta_expo         ! = K in:     B = A x Neff ** K
41
42real,dimension(nx,ny) :: neff   ! pression effective noeuds majeurs
43real,dimension(nx,ny) :: hwatmx ! hauteur d'eau staggered grille - afq
44real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq
45
46real :: coef_ile      ! coef frottement zones iles   (0.1)
47
48real,dimension(nx,ny) :: h_sedimmx   ! sediment thickness (or rebounded bedrock)
49real,dimension(nx,ny) :: h_sedimmy   ! sediment thickness (or rebounded bedrock)
50real :: seuil_sedim   ! drag reduced beyond this threshold
51real :: coef_sedim    ! drag reduction factor
52real :: neffmin = 1.e5 ! as in neffect, but this should be the exact same variable
53
54real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
55real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
56real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
57real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
58logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
59
60contains
61!-------------------------------------------------------------------------------
62
63!> SUBROUTINE: init_dragging
64!! Cette routine fait l'initialisation du dragging.
65!>
66subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
67
68implicit none
69
70real,dimension(nx,ny) :: h_sedim   ! sediment thickness (or rebounded bedrock)
71character(len=150) :: file_sedim   ! file for sediment thickness for HN or rebounded bsoc for Antar
72character(len=150) :: file_ncdf    ! fichier netcdf issue des fichiers .dat
73
74namelist/drag_param_beta_sedim/beta_slope,beta_expo,betamax,betamin,coef_ile,file_sedim,seuil_sedim,coef_sedim
75
76if (itracebug.eq.1)  call tracebug(' dragging param avec sedim')
77
78
79! formats pour les ecritures dans 42
80428 format(A)
81
82! lecture des parametres du run                      block drag neff
83!--------------------------------------------------------------------
84
85rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
86read(num_param,drag_param_beta_sedim)
87
88write(num_rep_42,428)'!___________________________________________________________' 
89write(num_rep_42,428) '&drag_param_beta_sedim          ! nom du bloc dragging param beta'
90write(num_rep_42,*)
91write(num_rep_42,*) 'beta_slope            = ', beta_slope
92write(num_rep_42,*) 'beta_expo             = ', beta_expo
93write(num_rep_42,*) 'betamax               = ', betamax
94write(num_rep_42,*) 'betamin               = ', betamin
95write(num_rep_42,*) 'file_sedim     = ', file_sedim
96write(num_rep_42,*) 'seuil_sedim    = ', seuil_sedim
97write(num_rep_42,*) 'coef_sedim     = ', coef_sedim
98write(num_rep_42,*)'/'                           
99write(num_rep_42,428) '! Slope & expo of beta = - slope x Neff ** expo'
100write(num_rep_42,428) '! Where h_sedim > seuil_sedim, beta*coef_sedim'
101
102!-------------------------------------------------------------------
103! masque stream
104      mstream_mx(:,:)=1
105      mstream_my(:,:)=1
106
107
108! coefficient permettant de modifier le basal drag.
109      drag_mx(:,:)=1.
110      drag_my(:,:)=1.
111     
112      betamax_2d(:,:) = betamax
113
114! modification of basal drag depends on where we have sediments
115! for Antarctica, we assume that what is present-day below sea level has sediment
116      file_sedim=trim(dirnameinp)//trim(file_sedim)
117      call lect_input(1,'h_sedim',1,h_sedim(:,:),file_sedim,file_ncdf)
118 
119      h_sedimmx(1,:)=h_sedim(1,:)
120      h_sedimmy(:,1)=h_sedim(:,1)
121      do i=2,nx
122         h_sedimmx(i,:)=(h_sedim(i-1,:)+h_sedim(i,:))/2.
123      enddo
124      do j=2,ny
125         h_sedimmy(:,j)=(h_sedim(:,j-1)+h_sedim(:,j))/2.
126      enddo
127 
128      toblim = 0.7e5 ! afq -- pour les iles, mais est-ce vraiment utile?
129
130      return
131      end subroutine init_dragging
132!________________________________________________________________________________
133
134!> SUBROUTINE: dragging
135!! Defini la localisation des streams et le frottement basal
136!>
137!-------------------------------------------------------------------------
138subroutine dragging   ! defini la localisation des streams et le frottement basal
139!$ USE OMP_LIB
140
141!         les iles n'ont pas de condition neff mais ont la condition toblim
142!         (ce bloc etait avant dans flottab)
143
144!$OMP PARALLEL
145!$OMP DO
146do j=2,ny
147   do i=2,nx
148      ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim)
149      ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim)
150   end do
151end do
152!$OMP END DO
153
154
155! le gzmx initial a ete calculé dans flottab pour les points a cheval sur la grounding line
156
157
158!$OMP WORKSHARE
159fleuvemx(:,:)=.false.
160fleuvemy(:,:)=.false.
161cote(:,:)=.false.
162gzmx(:,:)=.true.
163gzmy(:,:)=.true.
164flgzmx(:,:)=.false.
165flgzmy(:,:)=.false.
166!$OMP END WORKSHARE
167
168! detection des cotes
169!$OMP DO
170do  j=2,ny-1 
171   do i=2,nx-1 
172      if ((.not.flot(i,j)).and. & 
173      ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then
174         cote(i,j)=.true.
175      endif
176   end do
177end do
178!$OMP END DO
179
180! calcul de neff (pression effective sur noeuds majeurs)
181if (sum(neffmx(:,:)).le.0.) neffmx(:,:) =1.e8
182if (sum(neffmy(:,:)).le.0.) neffmy(:,:) =1.e8
183
184!$OMP DO
185do j=1,ny-1
186   do i=1,nx-1
187        neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
188   enddo
189enddo
190!$OMP END DO
191!aurel, for the borders:
192!$OMP WORKSHARE
193neff(:,ny)=neffmin
194neff(nx,:)=neffmin
195! calcul de hwat (staggered)
196hwatmx(:,:)=0.
197hwatmy(:,:)=0.
198!$OMP END WORKSHARE
199!$OMP DO
200do i=2,nx
201   hwatmx(i,:)=(hwater(i-1,:)+hwater(i,:))/2.
202enddo
203!$OMP END DO
204!$OMP DO
205do j=2,ny
206   hwatmy(:,j)=(hwater(:,j-1)+hwater(:,j))/2.
207enddo
208!$OMP END DO
209
210!$OMP WORKSHARE
211
212where (h_sedimmx(:,:).le.seuil_sedim)
213  betamx(:,:)= beta_slope*(neffmx(:,:)**beta_expo)
214elsewhere
215  betamx(:,:)= beta_slope*(neffmx(:,:)**beta_expo) * coef_sedim
216endwhere
217
218where (h_sedimmy(:,:).le.seuil_sedim)
219  betamy(:,:)= beta_slope*(neffmy(:,:)**beta_expo)
220elsewhere
221  betamy(:,:)= beta_slope*(neffmy(:,:)**beta_expo) * coef_sedim
222endwhere
223!where (h_sedimmx(:,:).gt.seuil_sedim) betamx(:,:) = beta_slope*(neffmin**beta_expo)/neffmin
224!where (h_sedimmy(:,:).gt.seuil_sedim) betamy(:,:) = beta_slope*(neffmin**beta_expo)/neffmin
225
226where (ilemx(:,:)) betamx(:,:) = betamx(:,:) * coef_ile
227where (ilemy(:,:)) betamy(:,:) = betamy(:,:) * coef_ile
228
229betamx(:,:)=max(betamx(:,:),betamin)
230betamy(:,:)=max(betamy(:,:),betamin)
231betamx(:,:)=min(betamx(:,:),betamax)
232betamy(:,:)=min(betamy(:,:),betamax)
233
234where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax
235where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax
236
237!$OMP END WORKSHARE
238
239
240! calcul de gzmx
241
242! points flottants : flgzmx mais pas gzmx
243!$OMP DO
244do j=2,ny
245   do i=2,nx
246      if (flot(i,j).and.(flot(i-1,j))) then
247         flotmx(i,j)=.true.
248         flgzmx(i,j)=.true.
249         gzmx(i,j)=.false.
250         betamx(i,j)=0.
251      end if
252      if (flot(i,j).and.(flot(i,j-1))) then
253         flotmy(i,j)=.true.
254         flgzmy(i,j)=.true.
255         gzmy(i,j)=.false.
256         betamy(i,j)=0.
257      end if
258   end do
259end do
260!$OMP END DO
261
262!--------- autres criteres
263
264! points poses
265!$OMP WORKSHARE
266gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.betamax)   !  Pas de calcul pour les points
267gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.betamax)   !  au fort beta
268
269flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
270flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
271
272fleuvemx(:,:)=gzmx(:,:)
273fleuvemy(:,:)=gzmy(:,:)
274!$OMP END WORKSHARE
275
276!$OMP DO
277do j=2,ny-1
278   do i=2,nx-1
279      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
280          + (betamy(i,j)+betamy(i,j+1)))*0.25
281   end do
282end do
283!$OMP END DO
284!$OMP END PARALLEL
285
286return
287end subroutine dragging
288
289end module dragging_param_beta_sedim
290
Note: See TracBrowser for help on using the repository browser.