source: branches/iLoveclim/SOURCES/spinup_mod.f90 @ 187

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

Grisli-iloveclim branch merged to trunk at revision 185
File size: 22.7 KB
Line 
1!> \file spinup_mod.f90
2!! This module allows various spinup tasks mostly based on balance velocities.
3!<
4
5!> \namespace spinup_vitbil
6!! This module allows various spinup tasks mostly based on balance velocities.
7!! @note Could be transformed to use observed surface velocities
8!! @note must be used as a dragging module 
9!! @note in module "choix" other dragging module must be switched off
10!! \author ...
11!! \date ...
12!! @note Used modules:
13!! @note  - use module3D_phy
14!! @note  - use param_phy_mod
15!! @note  - use deformation_mod_2lois
16!<
17
18
19module spinup_vitbil
20
21  use module3D_phy
22  use param_phy_mod
23  use deformation_mod_2lois
24  use interface_input
25  use io_netcdf_grisli
26 
27  implicit none
28
29  real, dimension(nx,ny) :: Vcol_x           !< vertically averaged velocity x direction (balance)
30  real, dimension(nx,ny) :: Vcol_y           !< vertically averaged velocity y direction (balance)
31  real, dimension(nx,ny) :: Ux_deformation   !< vertically averaged deformation velocity x direction
32  real, dimension(nx,ny) :: Uy_deformation   !< vertically averaged deformation velocity y direction
33  real, dimension(nx,ny) :: coef_defmx       !< rescaling coefficient of Sa_mx and s2a_mx
34  real, dimension(nx,ny) :: coef_defmy       !< rescaling coefficient of Sa_my and s2a_my
35  real, dimension(nx,ny) :: init_coefmx      !< rescalling coefficient before limitation
36  real, dimension(nx,ny) :: init_coefmy      !< rescalling coefficient before limitation
37
38! compatibilites avec remplimat
39  logical :: corr_def = .False.               !< for deformation correction (compatibility)
40
41  real, dimension(nx,ny) :: Vsl_x             !<
42  real, dimension(nx,ny) :: Vsl_y             !<
43 
44  integer :: type_vitbil ! defini le type de fichier a lire : lect_vitbil ou lect_vitbil_Lebrocq
45
46contains
47  !--------------------------------------------------------------------------------
48
49  !> SUBROUTINE: init_spinup
50  !! Initialize spinup
51  !<
52  subroutine init_spinup
53  namelist/spinup/ispinup,type_vitbil
54
55    ! put here which type of spinup
56    ! ispinup = 0                          ! run standard
57    ! ispinup = 1                          ! ispinup = 1 -> saute icethick3, diagnoshelf,
58    ! diffusiv, isostasie   pour equilibre temperature
59    ! en prenant les vitesses calculees aux premiers
60    ! pas de temps 
61    ! ispinup = 2                          ! ispinup = 2 -> fait la conservation
62    !                                      ! de la masse avec les vitesses de bilan
63    ! ispinup = 3                          ! fait le calcul des temperatures
64    !                                      ! avec les vitesses de bilan
65 
66
67
68
69! lecture des parametres
70
71rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
72read(num_param,spinup)
73
74write(num_rep_42,*)'!___________________________________________________________'
75write(num_rep_42,*)'!        spinup          module spinup_vitbil               '
76write(num_rep_42,spinup)                   ! pour ecrire les parametres lus
77write(num_rep_42,*)
78write(num_rep_42,*)
79write(num_rep_42,*)'! ispinup = 0     run standard voir no_spinup'
80write(num_rep_42,*)'! ispinup = 1     equilibre temperature avec vitesses grisli voir no_spinup'
81write(num_rep_42,*)
82write(num_rep_42,*)'! ispinup >1      use spinup_vitbil'
83write(num_rep_42,*)'! ispinup = 2     conservation de la masse avec vitesses bilan '
84write(num_rep_42,*)'! ispinup = 3     equilibre temperature avec vitesses bilan'
85write(num_rep_42,*)'! type_vitbil type de lecture des vitesses de bilan 1: lect_vitbil, 2 : lect_vitbil_Lebrocq'
86write(num_rep_42,*)
87
88 
89if (ispinup.eq.0) then
90   write(6,*)' ispinup = 0 et ispinup = 1 doivent etre appeles par no_spinup '
91   write(6,*) 'et il faut rajouter un dragging'
92   stop
93endif
94
95select case (type_vitbil)
96case (1)
97!cdc Methode Cat lecture fichier selon x et y sur grille stag
98  call lect_vitbil
99case(2)
100!cdc Methode Tof lecture fichier Lebrocq vitesse et direction
101  call lect_vitbil_Lebrocq
102case default
103  print*,'type_vitbil valeur invalide dans spinup_vitbil'
104  stop
105end select
106
107    if (itracebug.eq.1)  call tracebug(' fin routine init_spinup de spinup_vitbil')
108  end subroutine init_spinup
109
110  subroutine lect_vitbil
111
112    character(len=100) :: balance_Ux_file        !  balance velocity file Ux
113    character(len=100) :: balance_Uy_file        !  balance velocity file Uy
114          real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
115         
116    namelist/vitbil_upwind/balance_Ux_file, balance_Uy_file 
117
118    if (itracebug.eq.1)  call tracebug(' Subroutine lect_vitbil') 
119
120    ! lecture des parametres du run                   
121    !--------------------------------------------------------------------
122
123    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
124428 format(A)
125    read(num_param,vitbil_upwind)
126
127    write(num_rep_42,428) '!___________________________________________________________' 
128    write(num_rep_42,428) '!read balance velocities on staggered grid'
129    write(num_rep_42,vitbil_upwind)                   
130    write(num_rep_42,428) '! balance_Ux_file : nom du fichier qui contient les vit. bilan Ux'
131    write(num_rep_42,428) '! balance_Uy_file : nom du fichier qui contient les vit. bilan Uy'
132    write(num_rep_42,*)
133
134    ! read balance velocities on staggered nodes
135   
136    balance_Ux_file = trim(dirnameinp)//trim(balance_Ux_file)
137    balance_Uy_file = trim(dirnameinp)//trim(balance_Uy_file)
138   
139    call Read_Ncdf_var('Vcol_x',balance_Ux_file,tab)
140    Vcol_x(:,:)=tab(:,:)
141    call Read_Ncdf_var('Vcol_y',balance_Uy_file,tab)
142    Vcol_y(:,:)=tab(:,:)
143
144!    call lect_input(2,'Vcol',1,Vcol_x,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_x
145!    call lect_input(2,'Vcol',2,Vcol_y,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_y
146    !call lect_datfile(nx,ny,Vcol_x,1,balance_vel_file)   ! read Vcol_x   
147    !call lect_datfile(nx,ny,Vcol_y,2,balance_vel_file)   ! read Vcol_y
148
149    debug_3D(:,:,59)=Vcol_x(:,:)
150    debug_3D(:,:,60)=Vcol_y(:,:)
151
152  end subroutine lect_vitbil
153 
154  subroutine lect_vitbil_Lebrocq
155
156
157    character(len=100) :: vit_balance_file        !  balance velocity file
158    character(len=100) :: flowdir_balance_file    !  balance flow direction file
159    real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
160    real, dimension(nx,ny) :: V_Lebrocq           !< vertically averaged velocity (balance)
161    real, dimension(nx,ny) :: flowdir_Lebrocq     !< flow direction (degree)
162    real, dimension(nx,ny) :: Vx_Lebrocq, Vy_Lebrocq ! vitesses selon x et y
163   
164    namelist/vitbil_Lebrocq/vit_balance_file, flowdir_balance_file 
165
166    if (itracebug.eq.1)  call tracebug(' Subroutine lect_vitbil_Lebrocq') 
167
168    ! lecture des parametres du run                   
169    !--------------------------------------------------------------------
170
171    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
172428 format(A)
173    read(num_param,vitbil_Lebrocq)
174
175    write(num_rep_42,428) '!___________________________________________________________' 
176    write(num_rep_42,428) '!read balance velocities from Lebrocq code'
177    write(num_rep_42,vitbil_Lebrocq)                   
178    write(num_rep_42,428) '! vit_balance_file : balance velocity file'
179    write(num_rep_42,428) '! flowdir_balance_file : balance flow direction file'
180    write(num_rep_42,*)
181
182    ! read balance velocities on staggered nodes
183   
184    vit_balance_file = trim(dirnameinp)//trim(vit_balance_file)
185    flowdir_balance_file = trim(dirnameinp)//trim(flowdir_balance_file)
186   
187    call Read_Ncdf_var('z',vit_balance_file,tab)
188    V_Lebrocq(:,:)=tab(:,:)
189    call Read_Ncdf_var('z',flowdir_balance_file,tab)
190    flowdir_Lebrocq(:,:)=tab(:,:)*PI/180.
191
192                ! calcul des vitesses selon x et selon y :
193                where (V_Lebrocq(:,:).ge.0.)
194                        Vx_Lebrocq(:,:) = V_Lebrocq(:,:)*sin(flowdir_Lebrocq(:,:))
195                        Vy_Lebrocq(:,:) = V_Lebrocq(:,:)*cos(flowdir_Lebrocq(:,:))
196                elsewhere
197                        Vx_Lebrocq(:,:) = 0.
198                        Vy_Lebrocq(:,:) = 0.
199                endwhere
200               
201                ! calcul des vitesses sur les points staggered
202               
203                do j=2,ny
204                        do i=2,nx
205                                Vcol_x(i,j) = (Vx_Lebrocq(i-1,j) + Vx_Lebrocq(i,j))/2.
206                                Vcol_y(i,j) = (Vy_Lebrocq(i,j-1) + Vy_Lebrocq(i,j))/2.
207                        enddo
208                enddo
209
210  end subroutine lect_vitbil_Lebrocq 
211  !____________________________________________________________________________________
212
213  !> SUBROUTINE: init_dragging
214  !! @ Cette routine fait l'initialisation du dragging.
215  !<
216  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
217
218    implicit none
219
220    mstream_mx(:,:)=0             ! pas de dragging a l'initialisation
221    mstream_my(:,:)=0
222
223    ! coefficient permettant de modifier le basal drag.
224    drag_mx(:,:)=1.
225    drag_my(:,:)=1.
226    gzmx(:,:)=.false.
227    gzmy(:,:)=.false.
228    flgzmx(:,:)=flotmx(:,:)
229    flgzmy(:,:)=flotmy(:,:)
230
231  end subroutine init_dragging
232
233
234  subroutine dragging
235
236  end subroutine dragging
237  !----------------------------------------------------------------------
238
239
240  subroutine force_balance_vel
241
242    if (itracebug.eq.1)  call tracebug(' Subroutine force_balance_vel')
243
244    Uxbar(:,:)=Vcol_x(:,:)
245    Uybar(:,:)=Vcol_y(:,:)
246    debug_3D(:,:,59)=Uxbar(:,:)
247    debug_3D(:,:,60)=Uybar(:,:)
248  end subroutine force_balance_vel
249
250
251
252  !> SUBROUTINE: calc_coef_vitbil
253  !! @ Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way.
254  !! @note Must be called after flowlaw and before velocities 
255  !<
256  subroutine calc_coef_vitbil
257
258    ddbx(:,:)=0.
259    ddby(:,:)=0.
260    gzmx(:,:)=.true.
261    gzmy(:,:)=.true.
262    flgzmx(:,:)=.false.
263    flgzmy(:,:)=.false.
264    fleuvemx(:,:)=.false.
265    fleuvemy(:,:)=.false.
266   
267    do j=2,ny
268      do i=2,nx
269          if (flot(i,j).and.(flot(i-1,j))) then
270            flotmx(i,j)=.true.
271            flgzmx(i,j)=.true.
272            gzmx(i,j)=.false.
273          end if
274          if (flot(i,j).and.(flot(i,j-1))) then
275            flotmy(i,j)=.true.
276            flgzmy(i,j)=.true.
277            gzmy(i,j)=.false.
278          end if
279      end do
280    end do   
281   
282    where (coefmxbmelt(:,:).le.0.)
283      gzmx(:,:)=.false.
284    endwhere
285    where (coefmybmelt(:,:).le.0.)
286      gzmy(:,:)=.false.
287    endwhere
288    flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
289    flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
290    fleuvemx(:,:)=gzmx(:,:)
291    fleuvemy(:,:)=gzmy(:,:)
292   
293    if (itracebug.eq.1)  call tracebug(' Subroutine  calc_coef_vitbil')
294    call slope_surf
295
296    where (abs(sdx(:,:)).lt.1.e-6)
297       sdx(:,:)=-sign(1.e-6,Vcol_x(:,:))
298    end where
299    where (abs(sdy(:,:)).lt.1.e-6)
300       sdy(:,:)=-sign(1.e-6,Vcol_y(:,:))
301    end where
302
303
304    call calc_ubar_def
305
306    if (itracebug.eq.1)  call tracebug(' apres calc_ubar_def')
307
308    !---------------compute coef_defmx -----------------------------------
309
310    do j=1,ny
311       do i=1,nx
312
313          flotx:      if (flotmx(i,j)) then
314
315             coef_defmx(i,j) = 1.
316             Uxbar(i,j)      = Vcol_x(i,j)
317! dmr below is a cast from a real*4 to logical*4
318! dmr cannot be implicit in gfortran
319! dmr              flgzmx(i,j)     = Vcol_x(i,j)
320             flgzmx(i,j)     = (nint(Vcol_x(i,j)).ne.0)
321             uxdef(i,j)      = 0.
322             Ubx(i,j)        = Vcol_x(i,j)
323
324          else
325             coldx:      if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then   !base froide
326                !meme test que dans sliding Bindshadler
327                ! sans doute trop fort
328                ddbx(i,j) = 0.
329                Ubx(i,j)  = 0.
330
331                if (abs(Ux_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
332                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
333
334                else
335                   coef_defmx(i,j)=1.
336                end if
337             else                                                      ! base temperee
338                if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation
339                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
340                   ddbx(i,j) = 0
341                    Ubx(i,j)  = 0.
342                else
343                   coef_defmx(i,j)=1.                                  ! no rescaling of deformation
344                   ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)    ! pb si directions opposees
345                   Ubx(i,j) = Vcol_x(i,j)-Ux_deformation(i,j)
346                end if
347             end if coldx
348
349          end if flotx
350
351       end do
352    end do
353
354
355    do j=1,ny
356       do i=1,nx
357
358          floty:      if (flotmy(i,j)) then
359             coef_defmy(i,j) = 1.
360             Uybar(i,j)      = Vcol_y(i,j)
361! dmr below is a cast from a real*4 to logical*4
362! dmr cannot be implicit in gfortran
363! dmr             flgzmy(i,j)     = Vcol_y(i,j)
364             flgzmy(i,j)     = (nint(Vcol_y(i,j)).ne.0)
365
366             uydef(i,j)      = 0.
367             Uby(i,j)        = Vcol_y(i,j)
368
369 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test0')
370          else
371             coldy:      if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then   !base froide
372                !meme test que dans sliding Bindshadler
373                ddby(i,j) = 0.
374                Uby(i,j)  = 0.
375
376                if (abs(Uy_deformation(i,j)).gt.0.01) then             ! calcule par calc_ubar_def
377                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
378 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test1')
379                else
380 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test2')
381                   coef_defmy(i,j)=1.
382                end if
383             else                                                      ! base temperee
384                if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! que deformation
385                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
386                   ddby(i,j) = 0.
387                   Uby(i,j)  = 0.
388 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test3')
389                else
390                   coef_defmy(i,j)=1.                                  ! pas de rescaling de la deformation
391                   ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)        ! pb si directions opposees
392                   Uby(i,j) =  Vcol_y(i,j)-Uy_deformation(i,j)
393if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144))  call tracebug(' test4')
394
395                end if
396             end if coldy
397
398          end if floty
399
400       end do
401    end do
402 if (itracebug.eq.1)  call tracebug(' apres test floty')
403
404    ! remarque coef peut être <1  a cause de la pente locale, mais il joue toujours son role
405    debug_3D(:,:,61)=coef_defmx(:,:)
406    debug_3D(:,:,62)=coef_defmy(:,:)
407
408
409    calib: do iglen=n1poly,n2poly
410       do k=1,nz
411          SA_mx(:,:,k,iglen)  = SA_mx(:,:,k,iglen)*coef_defmx(:,:)
412          S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:)
413          SA_my(:,:,k,iglen)  = SA_my(:,:,k,iglen)*coef_defmy(:,:)
414          S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:)
415       end do
416    end do calib
417
418  end subroutine calc_coef_vitbil
419
420  !> SUBROUTINE: limit_coef_vitbil
421  !! Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way.
422  !! @note The difference with calc_coef_vitbil is that the coefficient is limited to the range 0.5-2
423  !! @note  - if > 2, sliding is assumed.
424  !! @note  - if < 0.5, it means that ice is actually colder than computed (to be implemented)
425  !! @note  Must be called after flowlaw and before velocities
426  !<
427
428  subroutine limit_coef_vitbil
429
430    call slope_surf
431
432    where (abs(sdx(:,:)).lt.1.e-6)
433       sdx(:,:)=-sign(1.e-6,Vcol_x(:,:))
434    end where
435    where (abs(sdy(:,:)).lt.1.e-6)
436       sdy(:,:)=-sign(1.e-6,Vcol_y(:,:))
437    end where
438
439
440    call calc_ubar_def
441
442    if (itracebug.eq.1)  call tracebug(' Entree dans routine calc_coef_vitbil') 
443    !---------------compute coef_defmx -----------------------------------
444
445    do j=1,ny
446       do i=1,nx
447
448          flotx:      if (flotmx(i,j)) then
449             coef_defmx(i,j)=1.
450             init_coefmx(i,j)= coef_defmx(i,j)
451          else
452             coldx:      if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then   !base froide
453                !meme test que dans sliding Bindshadler
454                ! sans doute trop fort
455                ddbx(i,j)=0.
456                if (abs(Ux_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
457                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
458                   init_coefmx(i,j)= coef_defmx(i,j)
459
460                   if  (coef_defmx(i,j).gt.1.) then
461                      coef_defmx(i,j)=1.
462                      ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)
463                   end if
464                else
465                   coef_defmx(i,j)=1.
466                   init_coefmx(i,j)= coef_defmx(i,j)
467                end if
468             else                                                      ! base temperee
469                if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation
470                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
471                   init_coefmx(i,j)= coef_defmx(i,j)
472                   ddbx(i,j) = 0
473                else
474                   coef_defmx(i,j)=1.                                  ! no rescaling of deformation
475                   init_coefmx(i,j)= coef_defmx(i,j)
476                   ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)    ! pb si directions opposees
477                end if
478             end if coldx
479
480          end if flotx
481
482       end do
483    end do
484
485
486    do j=1,ny
487       do i=1,nx
488
489          floty:      if (flotmy(i,j)) then
490             coef_defmy(i,j)=1.
491             init_coefmy(i,j)= coef_defmy(i,j)
492          else
493             coldy:      if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then   !base froide
494                !meme test que dans sliding Bindshadler
495                ! sans doute trop fort
496                ddby(i,j)=0.
497                if (abs(Uy_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
498                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
499                   init_coefmy(i,j)= coef_defmy(i,j)
500
501                   if  (coef_defmy(i,j).gt.1.) then
502                      coef_defmy(i,j)=1.
503                      ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)
504                   end if
505                else
506                   coef_defmy(i,j)=1.
507                   init_coefmy(i,j)= coef_defmy(i,j)
508                end if
509             else                                                      ! base temperee
510                if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! only deformation
511                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
512                   init_coefmy(i,j)= coef_defmy(i,j)
513                   ddby(i,j) = 0
514                else
515                   coef_defmy(i,j)=1.                                  ! no rescaling of deformation
516                   init_coefmy(i,j)= coef_defmy(i,j)
517                   ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)    ! pb si directions opposees
518                end if
519             end if coldy
520
521          end if floty
522
523       end do
524    end do
525
526
527    ! remarque coef peut être <1  a cause de la pente locale, mais il joue toujours son role
528    debug_3D(:,:,73)=init_coefmx(:,:)
529    debug_3D(:,:,74)=init_coefmy(:,:)
530
531    debug_3D(:,:,61)=coef_defmx(:,:)
532    debug_3D(:,:,62)=coef_defmy(:,:)
533
534
535    calib: do iglen=n1poly,n2poly
536       do k=1,nz
537          SA_mx(:,:,k,iglen)  = SA_mx(:,:,k,iglen)*coef_defmx(:,:)
538          S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:)
539          SA_my(:,:,k,iglen)  = SA_my(:,:,k,iglen)*coef_defmy(:,:)
540          S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:)
541       end do
542    end do calib
543
544    ! call ajust_ghf  ATTENTION NE PAS ACTIVER SAUF TESTS
545
546    !   debug_3D(:,:,73)= init_coefmx(:,:)
547    !   debug_3D(:,:,74)= init_coefmy(:,:)
548
549
550  end subroutine limit_coef_vitbil
551
552
553  !> SUBROUTINE: calc_ubar_def
554  !! Calculate velocity due to deformation before calibration
555  !<
556
557  subroutine calc_ubar_def   ! calculate velocity due to deformation before calibration
558
559    implicit none              ! extrait de diffusiv pour calculer la partie deformation
560
561    real :: glenexp
562    real :: inv_4dx ! inverse de dx pour eviter les divisions =1/(4*dx)
563    real :: inv_4dy ! inverse de dy pour eviter les divisions =1/(4*dy)
564
565    inv_4dx=1./(4*dx)
566    inv_4dy=1./(4*dy)
567
568
569    do iglen=n1poly,n2poly
570       glenexp=max(0.,(glen(iglen)-1.)/2.)
571
572       do j=1,ny
573          do i=1,nx
574
575             if (.not.flotmy(i,j)) then
576                ddy(i,j,iglen)=((slope2my(i,j)**  &  ! SLOPE2my calcule dans slope_surf
577                     glenexp)*(rog)**glen(iglen)) &
578                     *Hmy(i,j)**(glen(iglen)+1)
579
580             endif
581             if (.not.flotmx(i,j)) then                           
582                ddx(i,j,iglen)=((slope2mx(i,j)**  &  ! slope2mx calcule en debut de routine
583                     glenexp)*(rog)**glen(iglen)) &
584                     *Hmx(i,j)**(glen(iglen)+1)
585             endif
586
587          end do
588       end do
589    end do
590
591    do j=2,ny-1
592       do i=2,nx-1
593          ux_deformation(i,j)=0.
594          Uy_deformation(i,j)=0.
595          do iglen=n1poly,n2poly
596             ux_deformation(i,j)=ux_deformation(i,j)+ddx(i,j,iglen)*s2a_mx(i,j,1,iglen)
597             Uy_deformation(i,j)=Uy_deformation(i,j)+ddy(i,j,iglen)*s2a_my(i,j,1,iglen)
598          end do
599       end do
600    end do
601
602    do j=2,ny-1
603       do i=2,nx-1
604          ux_deformation(i,j)=ux_deformation(i,j)*(-sdx(i,j))
605          Uy_deformation(i,j)=Uy_deformation(i,j)*(-sdy(i,j))
606       end do
607    end do
608    !debug_3D(:,:,63)=ux_deformation(:,:)
609    !debug_3D(:,:,64)=uy_deformation(:,:)
610
611    if (itracebug.eq.1)  call tracebug(' fin de calc_ubar_def ')
612
613
614  end subroutine calc_ubar_def
615
616  !> SUBROUTINE: ajust_ghf
617  !! Ajuste le flux geothermique pour avoir une temperature coherente
618  !! avec les vitesses de bilan
619  !<
620
621  subroutine ajust_ghf
622
623    implicit none
624    real,dimension(nx,ny) :: coefdef_maj     !< coefficient deformation
625    real :: increment_ghf
626    real :: ghf_min
627    increment_ghf=0.00001                    !< exprime en  W/m2   
628    ghf_min=0.025                                           
629
630    do j=2,ny-1
631       do i=2,nx-1
632          if (.not. flot(i,j)) then
633             coefdef_maj(i,j)= (init_coefmx(i,j)+init_coefmx(i+1,j))+ &
634                  (init_coefmy(i,j)+init_coefmy(i,j+1))
635             coefdef_maj(i,j)=0.25*coefdef_maj(i,j)
636
637             ! un coef trop grand indique eventuellement une base trop froide
638             if ((coefdef_maj(i,j).gt.1.5).and.(ibase(i,j).eq.1)) then   ! en base froide
639                ghf(i,j)=ghf(i,j)-SECYEAR*increment_ghf   ! attention ghf est negatif
640
641                ! un coef trop petit indique une base trop chaude
642             else if ((coefdef_maj(i,j).lt.0.7).and.(ghf(i,j).lt.-secyear*ghf_min)) then
643                ghf(i,j)=ghf(i,j)+SECYEAR*increment_ghf   ! attention ghf est negatif
644             end if
645          end if
646       end do
647    end do
648
649    debug_3D(:,:,75)=(ghf0(:,:)-ghf(:,:))*1000./secyear
650
651  end subroutine ajust_ghf
652
653
654
655end module spinup_vitbil
656
657
658
659
Note: See TracBrowser for help on using the repository browser.