source: branches/GRISLIv3/SOURCES/Draggings_modules/dragging_param_beta_sedim_mod.f90 @ 470

Last change on this file since 470 was 470, checked in by aquiquet, 5 months ago

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