source: trunk/SOURCES/spinup_mod.f90 @ 25

Last change on this file since 25 was 22, checked in by roche, 9 years ago

Petites adaptations diverses du code pour compilation en gfortran. Ajout d un Makefile flexible a option pour choisir ifort ou gfortran.
File size: 18.9 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  implicit none
26
27  real, dimension(nx,ny) :: Vcol_x           !< vertically averaged velocity x direction (balance)
28  real, dimension(nx,ny) :: Vcol_y           !< vertically averaged velocity y direction (balance)
29  real, dimension(nx,ny) :: Ux_deformation   !< vertically averaged deformation velocity x direction
30  real, dimension(nx,ny) :: Uy_deformation   !< vertically averaged deformation velocity y direction
31  real, dimension(nx,ny) :: coef_defmx       !< rescaling coefficient of Sa_mx and s2a_mx
32  real, dimension(nx,ny) :: coef_defmy       !< rescaling coefficient of Sa_my and s2a_my
33  real, dimension(nx,ny) :: init_coefmx      !< rescalling coefficient before limitation
34  real, dimension(nx,ny) :: init_coefmy      !< rescalling coefficient before limitation
35
36! compatibilites avec remplimat
37  logical :: corr_def = .False.               !< for deformation correction (compatibility)
38
39  real, dimension(nx,ny) :: Vsl_x             !<
40  real, dimension(nx,ny) :: Vsl_y             !<
41
42contains
43  !--------------------------------------------------------------------------------
44
45  !> SUBROUTINE: init_spinup
46  !! Initialize spinup
47  !<
48  subroutine init_spinup
49  namelist/spinup/ispinup
50
51    ! put here which type of spinup
52    ! ispinup = 0                          ! run standard
53    ! ispinup = 1                          ! ispinup = 1 -> saute icethick3, diagnoshelf,
54    ! diffusiv, isostasie   pour equilibre temperature
55    ! en prenant les vitesses calculees aux premiers
56    ! pas de temps 
57    ! ispinup = 2                          ! ispinup = 2 -> fait la conservation
58    !                                      ! de la masse avec les vitesses de bilan
59    ! ispinup = 3                          ! fait le calcul des temperatures
60    !                                      ! avec les vitesses de bilan
61 
62
63
64
65! lecture des parametres
66
67rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
68read(num_param,spinup)
69
70write(num_rep_42,*)'!___________________________________________________________'
71write(num_rep_42,*)'!        spinup          module spinup_vitbil               '
72write(num_rep_42,spinup)                   ! pour ecrire les parametres lus
73write(num_rep_42,*)
74write(num_rep_42,*)
75write(num_rep_42,*)'! ispinup = 0     run standard voir no_spinup'
76write(num_rep_42,*)'! ispinup = 1     equilibre temperature avec vitesses grisli voir no_spinup'
77write(num_rep_42,*)
78write(num_rep_42,*)'! ispinup >1      use spinup_vitbil'
79write(num_rep_42,*)'! ispinup = 2     conservation de la masse avec vitesses bilan '
80write(num_rep_42,*)'! ispinup = 3     equilibre temperature avec vitesses bilan'
81write(num_rep_42,*)
82
83 
84if (ispinup.eq.0) then
85   write(6,*)' ispinup = 0 et ispinup = 1 doivent etre appeles par no_spinup '
86   write(6,*) 'et il faut rajouter un dragging'
87   stop
88endif
89
90
91    call lect_vitbil
92
93    if (itracebug.eq.1)  call tracebug(' fin routine init_spinup de spinup_vitbil')
94  end subroutine init_spinup
95
96  subroutine lect_vitbil
97
98
99    character(len=100) :: balance_Ux_file        !  balance velocity file Ux
100    character(len=100) :: balance_Uy_file        !  balance velocity file Uy
101
102    namelist/vitbil_upwind/balance_Ux_file, balance_Uy_file 
103
104    if (itracebug.eq.1)  call tracebug(' Subroutine lect_vitbil') 
105
106    ! lecture des parametres du run                   
107    !--------------------------------------------------------------------
108
109    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
110428 format(A)
111    read(num_param,vitbil_upwind)
112
113    write(num_rep_42,428) '!___________________________________________________________' 
114    write(num_rep_42,428) '!read balance velocities on staggered grid'
115    write(num_rep_42,vitbil_upwind)                   
116    write(num_rep_42,428) '! balance_Ux_file : nom du fichier qui contient les vit. bilan Ux'
117    write(num_rep_42,428) '! balance_Uy_file : nom du fichier qui contient les vit. bilan Uy'
118    write(num_rep_42,*)
119
120    ! read balance velocities on staggered nodes
121   
122    balance_Ux_file = trim(dirnameinp)//trim(balance_Ux_file)
123    balance_Uy_file = trim(dirnameinp)//trim(balance_Uy_file)
124
125    call lect_input(1,'Vcol_x',1,Vcol_x,balance_Ux_file,"")
126    call lect_input(1,'Vcol_y',1,Vcol_y,balance_Uy_file,"")
127
128!    call lect_input(2,'Vcol',1,Vcol_x,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_x
129!    call lect_input(2,'Vcol',2,Vcol_y,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_y
130    !call lect_datfile(nx,ny,Vcol_x,1,balance_vel_file)   ! read Vcol_x   
131    !call lect_datfile(nx,ny,Vcol_y,2,balance_vel_file)   ! read Vcol_y
132
133    debug_3D(:,:,59)=Vcol_x(:,:)
134    debug_3D(:,:,60)=Vcol_y(:,:)
135
136  end subroutine lect_vitbil
137  !____________________________________________________________________________________
138
139  !> SUBROUTINE: init_dragging
140  !! @ Cette routine fait l'initialisation du dragging.
141  !<
142  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
143
144    implicit none
145
146    mstream_mx(:,:)=0             ! pas de dragging a l'initialisation
147    mstream_my(:,:)=0
148
149    ! coefficient permettant de modifier le basal drag.
150    drag_mx(:,:)=1.
151    drag_my(:,:)=1.
152    gzmx(:,:)=.false.
153    gzmy(:,:)=.false.
154    flgzmx(:,:)=flotmx(:,:)
155    flgzmy(:,:)=flotmy(:,:)
156
157  end subroutine init_dragging
158
159
160  subroutine dragging
161
162  end subroutine dragging
163  !----------------------------------------------------------------------
164
165
166  subroutine force_balance_vel
167
168    if (itracebug.eq.1)  call tracebug(' Subroutine force_balance_vel')
169
170    Uxbar(:,:)=Vcol_x(:,:)
171    Uybar(:,:)=Vcol_y(:,:)
172    debug_3D(:,:,59)=Uxbar(:,:)
173    debug_3D(:,:,60)=Uybar(:,:)
174  end subroutine force_balance_vel
175
176
177
178  !> SUBROUTINE: calc_coef_vitbil
179  !! @ Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way.
180  !! @note Must be called after flowlaw and before velocities 
181  !<
182  subroutine calc_coef_vitbil
183
184    ddbx(:,:)=0.
185    ddby(:,:)=0.
186
187    if (itracebug.eq.1)  call tracebug(' Subroutine  calc_coef_vitbil')
188    call slope_surf
189
190    where (abs(sdx(:,:)).lt.1.e-6)
191       sdx(:,:)=-sign(1.e-6,Vcol_x(:,:))
192    end where
193    where (abs(sdy(:,:)).lt.1.e-6)
194       sdy(:,:)=-sign(1.e-6,Vcol_y(:,:))
195    end where
196
197
198    call calc_ubar_def
199
200    if (itracebug.eq.1)  call tracebug(' apres calc_ubar_def')
201
202    !---------------compute coef_defmx -----------------------------------
203
204    do j=1,ny
205       do i=1,nx
206
207          flotx:      if (flotmx(i,j)) then
208
209             coef_defmx(i,j) = 1.
210             Uxbar(i,j)      = Vcol_x(i,j)
211             flgzmx(i,j)     = transfer(Vcol_x(i,j),flgzmx(i,j))
212             uxdef(i,j)      = 0.
213             Ubx(i,j)        = Vcol_x(i,j)
214
215          else
216             coldx:      if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then   !base froide
217                !meme test que dans sliding Bindshadler
218                ! sans doute trop fort
219                ddbx(i,j) = 0.
220                Ubx(i,j)  = 0.
221
222                if (abs(Ux_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
223                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
224
225                else
226                   coef_defmx(i,j)=1.
227                end if
228             else                                                      ! base temperee
229                if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation
230                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
231                   ddbx(i,j) = 0
232                    Ubx(i,j)  = 0.
233                else
234                   coef_defmx(i,j)=1.                                  ! no rescaling of deformation
235                   ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)    ! pb si directions opposees
236                   Ubx(i,j) = Vcol_x(i,j)-Ux_deformation(i,j)
237                end if
238             end if coldx
239
240          end if flotx
241
242       end do
243    end do
244
245
246    do j=1,ny
247       do i=1,nx
248
249          floty:      if (flotmy(i,j)) then
250             coef_defmy(i,j) = 1.
251             Uybar(i,j)      = Vcol_y(i,j)
252             flgzmy(i,j)     = transfer(Vcol_y(i,j),flgzmy(i,j))
253             uydef(i,j)      = 0.
254             Uby(i,j)        = Vcol_y(i,j)
255
256 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test0')
257          else
258             coldy:      if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then   !base froide
259                !meme test que dans sliding Bindshadler
260                ddby(i,j) = 0.
261                Uby(i,j)  = 0.
262
263                if (abs(Uy_deformation(i,j)).gt.0.01) then             ! calcule par calc_ubar_def
264                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
265 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test1')
266                else
267 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test2')
268                   coef_defmy(i,j)=1.
269                end if
270             else                                                      ! base temperee
271                if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! que deformation
272                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
273                   ddby(i,j) = 0.
274                   Uby(i,j)  = 0.
275 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test3')
276                else
277                   coef_defmy(i,j)=1.                                  ! pas de rescaling de la deformation
278                   ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)        ! pb si directions opposees
279                   Uby(i,j) =  Vcol_y(i,j)-Uy_deformation(i,j)
280if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144))  call tracebug(' test4')
281
282                end if
283             end if coldy
284
285          end if floty
286
287       end do
288    end do
289 if (itracebug.eq.1)  call tracebug(' apres test floty')
290
291    ! remarque coef peut être <1  a cause de la pente locale, mais il joue toujours son role
292    debug_3D(:,:,61)=coef_defmx(:,:)
293    debug_3D(:,:,62)=coef_defmy(:,:)
294
295
296    calib: do iglen=n1poly,n2poly
297       do k=1,nz
298          SA_mx(:,:,k,iglen)  = SA_mx(:,:,k,iglen)*coef_defmx(:,:)
299          S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:)
300          SA_my(:,:,k,iglen)  = SA_my(:,:,k,iglen)*coef_defmy(:,:)
301          S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:)
302       end do
303    end do calib
304
305  end subroutine calc_coef_vitbil
306
307  !> SUBROUTINE: limit_coef_vitbil
308  !! Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way.
309  !! @note The difference with calc_coef_vitbil is that the coefficient is limited to the range 0.5-2
310  !! @note  - if > 2, sliding is assumed.
311  !! @note  - if < 0.5, it means that ice is actually colder than computed (to be implemented)
312  !! @note  Must be called after flowlaw and before velocities
313  !<
314
315  subroutine limit_coef_vitbil
316
317    call slope_surf
318
319    where (abs(sdx(:,:)).lt.1.e-6)
320       sdx(:,:)=-sign(1.e-6,Vcol_x(:,:))
321    end where
322    where (abs(sdy(:,:)).lt.1.e-6)
323       sdy(:,:)=-sign(1.e-6,Vcol_y(:,:))
324    end where
325
326
327    call calc_ubar_def
328
329    if (itracebug.eq.1)  call tracebug(' Entree dans routine calc_coef_vitbil') 
330    !---------------compute coef_defmx -----------------------------------
331
332    do j=1,ny
333       do i=1,nx
334
335          flotx:      if (flotmx(i,j)) then
336             coef_defmx(i,j)=1.
337             init_coefmx(i,j)= coef_defmx(i,j)
338          else
339             coldx:      if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then   !base froide
340                !meme test que dans sliding Bindshadler
341                ! sans doute trop fort
342                ddbx(i,j)=0.
343                if (abs(Ux_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
344                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
345                   init_coefmx(i,j)= coef_defmx(i,j)
346
347                   if  (coef_defmx(i,j).gt.1.) then
348                      coef_defmx(i,j)=1.
349                      ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)
350                   end if
351                else
352                   coef_defmx(i,j)=1.
353                   init_coefmx(i,j)= coef_defmx(i,j)
354                end if
355             else                                                      ! base temperee
356                if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation
357                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
358                   init_coefmx(i,j)= coef_defmx(i,j)
359                   ddbx(i,j) = 0
360                else
361                   coef_defmx(i,j)=1.                                  ! no rescaling of deformation
362                   init_coefmx(i,j)= coef_defmx(i,j)
363                   ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)    ! pb si directions opposees
364                end if
365             end if coldx
366
367          end if flotx
368
369       end do
370    end do
371
372
373    do j=1,ny
374       do i=1,nx
375
376          floty:      if (flotmy(i,j)) then
377             coef_defmy(i,j)=1.
378             init_coefmy(i,j)= coef_defmy(i,j)
379          else
380             coldy:      if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then   !base froide
381                !meme test que dans sliding Bindshadler
382                ! sans doute trop fort
383                ddby(i,j)=0.
384                if (abs(Uy_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
385                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
386                   init_coefmy(i,j)= coef_defmy(i,j)
387
388                   if  (coef_defmy(i,j).gt.1.) then
389                      coef_defmy(i,j)=1.
390                      ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)
391                   end if
392                else
393                   coef_defmy(i,j)=1.
394                   init_coefmy(i,j)= coef_defmy(i,j)
395                end if
396             else                                                      ! base temperee
397                if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! only deformation
398                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
399                   init_coefmy(i,j)= coef_defmy(i,j)
400                   ddby(i,j) = 0
401                else
402                   coef_defmy(i,j)=1.                                  ! no rescaling of deformation
403                   init_coefmy(i,j)= coef_defmy(i,j)
404                   ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)    ! pb si directions opposees
405                end if
406             end if coldy
407
408          end if floty
409
410       end do
411    end do
412
413
414    ! remarque coef peut être <1  a cause de la pente locale, mais il joue toujours son role
415    debug_3D(:,:,73)=init_coefmx(:,:)
416    debug_3D(:,:,74)=init_coefmy(:,:)
417
418    debug_3D(:,:,61)=coef_defmx(:,:)
419    debug_3D(:,:,62)=coef_defmy(:,:)
420
421
422    calib: do iglen=n1poly,n2poly
423       do k=1,nz
424          SA_mx(:,:,k,iglen)  = SA_mx(:,:,k,iglen)*coef_defmx(:,:)
425          S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:)
426          SA_my(:,:,k,iglen)  = SA_my(:,:,k,iglen)*coef_defmy(:,:)
427          S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:)
428       end do
429    end do calib
430
431    ! call ajust_ghf  ATTENTION NE PAS ACTIVER SAUF TESTS
432
433    !   debug_3D(:,:,73)= init_coefmx(:,:)
434    !   debug_3D(:,:,74)= init_coefmy(:,:)
435
436
437  end subroutine limit_coef_vitbil
438
439
440  !> SUBROUTINE: calc_ubar_def
441  !! Calculate velocity due to deformation before calibration
442  !<
443
444  subroutine calc_ubar_def   ! calculate velocity due to deformation before calibration
445
446    implicit none              ! extrait de diffusiv pour calculer la partie deformation
447
448    real :: glenexp
449    real :: inv_4dx ! inverse de dx pour eviter les divisions =1/(4*dx)
450    real :: inv_4dy ! inverse de dy pour eviter les divisions =1/(4*dy)
451
452    inv_4dx=1./(4*dx)
453    inv_4dy=1./(4*dy)
454
455
456    do iglen=n1poly,n2poly
457       glenexp=max(0.,(glen(iglen)-1.)/2.)
458
459       do j=1,ny
460          do i=1,nx
461
462             if (.not.flotmy(i,j)) then
463                ddy(i,j,iglen)=((slope2my(i,j)**  &  ! SLOPE2my calcule dans slope_surf
464                     glenexp)*(rog)**glen(iglen)) &
465                     *Hmy(i,j)**(glen(iglen)+1)
466
467             endif
468             if (.not.flotmx(i,j)) then                           
469                ddx(i,j,iglen)=((slope2mx(i,j)**  &  ! slope2mx calcule en debut de routine
470                     glenexp)*(rog)**glen(iglen)) &
471                     *Hmx(i,j)**(glen(iglen)+1)
472             endif
473
474          end do
475       end do
476    end do
477
478    do j=2,ny-1
479       do i=2,nx-1
480          ux_deformation(i,j)=0.
481          Uy_deformation(i,j)=0.
482          do iglen=n1poly,n2poly
483             ux_deformation(i,j)=ux_deformation(i,j)+ddx(i,j,iglen)*s2a_mx(i,j,1,iglen)
484             Uy_deformation(i,j)=Uy_deformation(i,j)+ddy(i,j,iglen)*s2a_my(i,j,1,iglen)
485          end do
486       end do
487    end do
488
489    do j=2,ny-1
490       do i=2,nx-1
491          ux_deformation(i,j)=ux_deformation(i,j)*(-sdx(i,j))
492          Uy_deformation(i,j)=Uy_deformation(i,j)*(-sdy(i,j))
493       end do
494    end do
495    debug_3D(:,:,63)=ux_deformation(:,:)
496    debug_3D(:,:,64)=uy_deformation(:,:)
497
498    if (itracebug.eq.1)  call tracebug(' fin de calc_ubar_def ')
499
500
501  end subroutine calc_ubar_def
502
503  !> SUBROUTINE: ajust_ghf
504  !! Ajuste le flux geothermique pour avoir une temperature coherente
505  !! avec les vitesses de bilan
506  !<
507
508  subroutine ajust_ghf
509
510    implicit none
511    real,dimension(nx,ny) :: corr_ghf        !< correction du gflux geothermique
512    real,dimension(nx,ny) :: coefdef_maj     !< coefficient deformation
513    real :: increment_ghf
514    real :: ghf_min
515    increment_ghf=0.00001                    !< exprime en  W/m2   
516    ghf_min=0.025                                           
517
518    do j=2,ny-1
519       do i=2,nx-1
520          if (.not. flot(i,j)) then
521             coefdef_maj(i,j)= (init_coefmx(i,j)+init_coefmx(i+1,j))+ &
522                  (init_coefmy(i,j)+init_coefmy(i,j+1))
523             coefdef_maj(i,j)=0.25*coefdef_maj(i,j)
524
525             ! un coef trop grand indique eventuellement une base trop froide
526             if ((coefdef_maj(i,j).gt.1.5).and.(ibase(i,j).eq.1)) then   ! en base froide
527                ghf(i,j)=ghf(i,j)-SECYEAR*increment_ghf   ! attention ghf est negatif
528
529                ! un coef trop petit indique une base trop chaude
530             else if ((coefdef_maj(i,j).lt.0.7).and.(ghf(i,j).lt.-secyear*ghf_min)) then
531                ghf(i,j)=ghf(i,j)+SECYEAR*increment_ghf   ! attention ghf est negatif
532             end if
533          end if
534       end do
535    end do
536
537    debug_3D(:,:,75)=(ghf0(:,:)-ghf(:,:))*1000./secyear
538
539  end subroutine ajust_ghf
540
541
542
543end module spinup_vitbil
544
545
546
547
Note: See TracBrowser for help on using the repository browser.