source: trunk/SOURCES/spinup_mod.f90 @ 29

Last change on this file since 29 was 29, checked in by dumas, 8 years ago

Mise a jour pour etre compatible avec le code couple iLOVECLIM
File size: 19.2 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! dmr below is a cast from a real*4 to logical*4
212! dmr cannot be implicit in gfortran
213! dmr              flgzmx(i,j)     = Vcol_x(i,j)
214             flgzmx(i,j)     = (nint(Vcol_x(i,j)).ne.0)
215             uxdef(i,j)      = 0.
216             Ubx(i,j)        = Vcol_x(i,j)
217
218          else
219             coldx:      if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then   !base froide
220                !meme test que dans sliding Bindshadler
221                ! sans doute trop fort
222                ddbx(i,j) = 0.
223                Ubx(i,j)  = 0.
224
225                if (abs(Ux_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
226                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
227
228                else
229                   coef_defmx(i,j)=1.
230                end if
231             else                                                      ! base temperee
232                if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation
233                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
234                   ddbx(i,j) = 0
235                    Ubx(i,j)  = 0.
236                else
237                   coef_defmx(i,j)=1.                                  ! no rescaling of deformation
238                   ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)    ! pb si directions opposees
239                   Ubx(i,j) = Vcol_x(i,j)-Ux_deformation(i,j)
240                end if
241             end if coldx
242
243          end if flotx
244
245       end do
246    end do
247
248
249    do j=1,ny
250       do i=1,nx
251
252          floty:      if (flotmy(i,j)) then
253             coef_defmy(i,j) = 1.
254             Uybar(i,j)      = Vcol_y(i,j)
255! dmr below is a cast from a real*4 to logical*4
256! dmr cannot be implicit in gfortran
257! dmr             flgzmy(i,j)     = Vcol_y(i,j)
258             flgzmy(i,j)     = (nint(Vcol_y(i,j)).ne.0)
259
260             uydef(i,j)      = 0.
261             Uby(i,j)        = Vcol_y(i,j)
262
263 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test0')
264          else
265             coldy:      if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then   !base froide
266                !meme test que dans sliding Bindshadler
267                ddby(i,j) = 0.
268                Uby(i,j)  = 0.
269
270                if (abs(Uy_deformation(i,j)).gt.0.01) then             ! calcule par calc_ubar_def
271                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
272 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test1')
273                else
274 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test2')
275                   coef_defmy(i,j)=1.
276                end if
277             else                                                      ! base temperee
278                if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! que deformation
279                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
280                   ddby(i,j) = 0.
281                   Uby(i,j)  = 0.
282 if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test3')
283                else
284                   coef_defmy(i,j)=1.                                  ! pas de rescaling de la deformation
285                   ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)        ! pb si directions opposees
286                   Uby(i,j) =  Vcol_y(i,j)-Uy_deformation(i,j)
287if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144))  call tracebug(' test4')
288
289                end if
290             end if coldy
291
292          end if floty
293
294       end do
295    end do
296 if (itracebug.eq.1)  call tracebug(' apres test floty')
297
298    ! remarque coef peut être <1  a cause de la pente locale, mais il joue toujours son role
299    debug_3D(:,:,61)=coef_defmx(:,:)
300    debug_3D(:,:,62)=coef_defmy(:,:)
301
302
303    calib: do iglen=n1poly,n2poly
304       do k=1,nz
305          SA_mx(:,:,k,iglen)  = SA_mx(:,:,k,iglen)*coef_defmx(:,:)
306          S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:)
307          SA_my(:,:,k,iglen)  = SA_my(:,:,k,iglen)*coef_defmy(:,:)
308          S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:)
309       end do
310    end do calib
311
312  end subroutine calc_coef_vitbil
313
314  !> SUBROUTINE: limit_coef_vitbil
315  !! Calibrate Sa and S2a to force velocity to be balance_velocity in a consistent way.
316  !! @note The difference with calc_coef_vitbil is that the coefficient is limited to the range 0.5-2
317  !! @note  - if > 2, sliding is assumed.
318  !! @note  - if < 0.5, it means that ice is actually colder than computed (to be implemented)
319  !! @note  Must be called after flowlaw and before velocities
320  !<
321
322  subroutine limit_coef_vitbil
323
324    call slope_surf
325
326    where (abs(sdx(:,:)).lt.1.e-6)
327       sdx(:,:)=-sign(1.e-6,Vcol_x(:,:))
328    end where
329    where (abs(sdy(:,:)).lt.1.e-6)
330       sdy(:,:)=-sign(1.e-6,Vcol_y(:,:))
331    end where
332
333
334    call calc_ubar_def
335
336    if (itracebug.eq.1)  call tracebug(' Entree dans routine calc_coef_vitbil') 
337    !---------------compute coef_defmx -----------------------------------
338
339    do j=1,ny
340       do i=1,nx
341
342          flotx:      if (flotmx(i,j)) then
343             coef_defmx(i,j)=1.
344             init_coefmx(i,j)= coef_defmx(i,j)
345          else
346             coldx:      if ((coefmxbmelt(i,j).le.0.).and.(.not.gzmx(i,j))) then   !base froide
347                !meme test que dans sliding Bindshadler
348                ! sans doute trop fort
349                ddbx(i,j)=0.
350                if (abs(Ux_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
351                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
352                   init_coefmx(i,j)= coef_defmx(i,j)
353
354                   if  (coef_defmx(i,j).gt.1.) then
355                      coef_defmx(i,j)=1.
356                      ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)
357                   end if
358                else
359                   coef_defmx(i,j)=1.
360                   init_coefmx(i,j)= coef_defmx(i,j)
361                end if
362             else                                                      ! base temperee
363                if (abs(Ux_deformation(i,j)).gt.abs(Vcol_x(i,j))) then ! only deformation
364                   coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
365                   init_coefmx(i,j)= coef_defmx(i,j)
366                   ddbx(i,j) = 0
367                else
368                   coef_defmx(i,j)=1.                                  ! no rescaling of deformation
369                   init_coefmx(i,j)= coef_defmx(i,j)
370                   ddbx(i,j)=-(Vcol_x(i,j)-Ux_deformation(i,j))/sdx(i,j)    ! pb si directions opposees
371                end if
372             end if coldx
373
374          end if flotx
375
376       end do
377    end do
378
379
380    do j=1,ny
381       do i=1,nx
382
383          floty:      if (flotmy(i,j)) then
384             coef_defmy(i,j)=1.
385             init_coefmy(i,j)= coef_defmy(i,j)
386          else
387             coldy:      if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then   !base froide
388                !meme test que dans sliding Bindshadler
389                ! sans doute trop fort
390                ddby(i,j)=0.
391                if (abs(Uy_deformation(i,j)).gt.0.01) then             ! calcule par  calc_ubar_def
392                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
393                   init_coefmy(i,j)= coef_defmy(i,j)
394
395                   if  (coef_defmy(i,j).gt.1.) then
396                      coef_defmy(i,j)=1.
397                      ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)
398                   end if
399                else
400                   coef_defmy(i,j)=1.
401                   init_coefmy(i,j)= coef_defmy(i,j)
402                end if
403             else                                                      ! base temperee
404                if (abs(Uy_deformation(i,j)).gt.abs(Vcol_y(i,j))) then ! only deformation
405                   coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
406                   init_coefmy(i,j)= coef_defmy(i,j)
407                   ddby(i,j) = 0
408                else
409                   coef_defmy(i,j)=1.                                  ! no rescaling of deformation
410                   init_coefmy(i,j)= coef_defmy(i,j)
411                   ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)    ! pb si directions opposees
412                end if
413             end if coldy
414
415          end if floty
416
417       end do
418    end do
419
420
421    ! remarque coef peut être <1  a cause de la pente locale, mais il joue toujours son role
422    debug_3D(:,:,73)=init_coefmx(:,:)
423    debug_3D(:,:,74)=init_coefmy(:,:)
424
425    debug_3D(:,:,61)=coef_defmx(:,:)
426    debug_3D(:,:,62)=coef_defmy(:,:)
427
428
429    calib: do iglen=n1poly,n2poly
430       do k=1,nz
431          SA_mx(:,:,k,iglen)  = SA_mx(:,:,k,iglen)*coef_defmx(:,:)
432          S2A_mx(:,:,k,iglen) = S2A_mx(:,:,k,iglen)*coef_defmx(:,:)
433          SA_my(:,:,k,iglen)  = SA_my(:,:,k,iglen)*coef_defmy(:,:)
434          S2A_my(:,:,k,iglen) = S2A_my(:,:,k,iglen)*coef_defmy(:,:)
435       end do
436    end do calib
437
438    ! call ajust_ghf  ATTENTION NE PAS ACTIVER SAUF TESTS
439
440    !   debug_3D(:,:,73)= init_coefmx(:,:)
441    !   debug_3D(:,:,74)= init_coefmy(:,:)
442
443
444  end subroutine limit_coef_vitbil
445
446
447  !> SUBROUTINE: calc_ubar_def
448  !! Calculate velocity due to deformation before calibration
449  !<
450
451  subroutine calc_ubar_def   ! calculate velocity due to deformation before calibration
452
453    implicit none              ! extrait de diffusiv pour calculer la partie deformation
454
455    real :: glenexp
456    real :: inv_4dx ! inverse de dx pour eviter les divisions =1/(4*dx)
457    real :: inv_4dy ! inverse de dy pour eviter les divisions =1/(4*dy)
458
459    inv_4dx=1./(4*dx)
460    inv_4dy=1./(4*dy)
461
462
463    do iglen=n1poly,n2poly
464       glenexp=max(0.,(glen(iglen)-1.)/2.)
465
466       do j=1,ny
467          do i=1,nx
468
469             if (.not.flotmy(i,j)) then
470                ddy(i,j,iglen)=((slope2my(i,j)**  &  ! SLOPE2my calcule dans slope_surf
471                     glenexp)*(rog)**glen(iglen)) &
472                     *Hmy(i,j)**(glen(iglen)+1)
473
474             endif
475             if (.not.flotmx(i,j)) then                           
476                ddx(i,j,iglen)=((slope2mx(i,j)**  &  ! slope2mx calcule en debut de routine
477                     glenexp)*(rog)**glen(iglen)) &
478                     *Hmx(i,j)**(glen(iglen)+1)
479             endif
480
481          end do
482       end do
483    end do
484
485    do j=2,ny-1
486       do i=2,nx-1
487          ux_deformation(i,j)=0.
488          Uy_deformation(i,j)=0.
489          do iglen=n1poly,n2poly
490             ux_deformation(i,j)=ux_deformation(i,j)+ddx(i,j,iglen)*s2a_mx(i,j,1,iglen)
491             Uy_deformation(i,j)=Uy_deformation(i,j)+ddy(i,j,iglen)*s2a_my(i,j,1,iglen)
492          end do
493       end do
494    end do
495
496    do j=2,ny-1
497       do i=2,nx-1
498          ux_deformation(i,j)=ux_deformation(i,j)*(-sdx(i,j))
499          Uy_deformation(i,j)=Uy_deformation(i,j)*(-sdy(i,j))
500       end do
501    end do
502    debug_3D(:,:,63)=ux_deformation(:,:)
503    debug_3D(:,:,64)=uy_deformation(:,:)
504
505    if (itracebug.eq.1)  call tracebug(' fin de calc_ubar_def ')
506
507
508  end subroutine calc_ubar_def
509
510  !> SUBROUTINE: ajust_ghf
511  !! Ajuste le flux geothermique pour avoir une temperature coherente
512  !! avec les vitesses de bilan
513  !<
514
515  subroutine ajust_ghf
516
517    implicit none
518    real,dimension(nx,ny) :: corr_ghf        !< correction du gflux geothermique
519    real,dimension(nx,ny) :: coefdef_maj     !< coefficient deformation
520    real :: increment_ghf
521    real :: ghf_min
522    increment_ghf=0.00001                    !< exprime en  W/m2   
523    ghf_min=0.025                                           
524
525    do j=2,ny-1
526       do i=2,nx-1
527          if (.not. flot(i,j)) then
528             coefdef_maj(i,j)= (init_coefmx(i,j)+init_coefmx(i+1,j))+ &
529                  (init_coefmy(i,j)+init_coefmy(i,j+1))
530             coefdef_maj(i,j)=0.25*coefdef_maj(i,j)
531
532             ! un coef trop grand indique eventuellement une base trop froide
533             if ((coefdef_maj(i,j).gt.1.5).and.(ibase(i,j).eq.1)) then   ! en base froide
534                ghf(i,j)=ghf(i,j)-SECYEAR*increment_ghf   ! attention ghf est negatif
535
536                ! un coef trop petit indique une base trop chaude
537             else if ((coefdef_maj(i,j).lt.0.7).and.(ghf(i,j).lt.-secyear*ghf_min)) then
538                ghf(i,j)=ghf(i,j)+SECYEAR*increment_ghf   ! attention ghf est negatif
539             end if
540          end if
541       end do
542    end do
543
544    debug_3D(:,:,75)=(ghf0(:,:)-ghf(:,:))*1000./secyear
545
546  end subroutine ajust_ghf
547
548
549
550end module spinup_vitbil
551
552
553
554
Note: See TracBrowser for help on using the repository browser.