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

Last change on this file since 504 was 504, checked in by aquiquet, 6 weeks ago

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