source: trunk/SOURCES/spinup_mod.f90 @ 334

Last change on this file since 334 was 206, checked in by dumas, 6 years ago

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