source: trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_mod.f90 @ 29

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

Mise a jour pour etre compatible avec le code couple iLOVECLIM

File size: 18.7 KB
Line 
1!> \file climat_GrIce2sea_mod.f90
2! forcage avec BM
3!! Module ou les variations temporelles des variables climatiques
4!! sont directement imposees
5!<
6
7!> \namespace  climat_Grice2sea_mod
8!! Module ou les variations temporelles des variables climatiques
9!! sont directement imposees
10!! \author Cat
11!! \date 31 oct
12!! @note Used modules:
13!! @note   - use module3D_phy
14!<
15module climat_Grice2sea_years_mod
16
17
18use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,dirnameinp
19!use lect_climref_Ice2sea
20use netcdf
21use io_netcdf_grisli
22implicit none
23
24
25real :: T_lapse_rate              !< pour la temperature
26
27
28! declarations specifiques year by year
29
30real,dimension(:),allocatable        :: time_snap         !> date des snapshots  indice : nb de nb_snap)
31real,dimension(:,:,:),allocatable    :: smb_snap          !> bilan de masse des snapshots indices : nx,ny,nb_snap
32real                                 :: ecart_snap        !> ecart temporel entre les snapshots
33real                                 :: time_depart_snaps !> temps du debut premier snapshot
34integer                              :: nb_snap           !> nombre de snapshots
35
36! declaration pour le bilan de masse
37real,dimension(nx,ny)                :: bm_0              !> bilan de masse initial equ. S_0ref de Tamsin
38real,dimension(nx,ny)                :: bm_time           !> bilan de masse interpole entre 2 snapshots (~ S_t^RCM)
39
40! calcul du gradient
41real,dimension(nx,ny)                :: bm_ref_t          !> bilan de masse de reference au temps t
42real,dimension(nx,ny,10)             :: bm_ref_10         !> tableau des bm_ref pour moyenne 10 ans
43real                                 :: time_prec         !> temps du precedent snapshot
44integer                              :: icum              !> pour le test stockage dans bm_ref_10
45integer                              :: i_moy             !> nombre de snapshots stockes
46
47real,dimension(nx,ny)                :: grad_bm           !> gradient du bilan de masse
48real,dimension(nx,ny)                :: S_ref             !> surface de reference
49
50real                                 :: grad_N_smb_neg    !> SMB < 0 North
51real                                 :: grad_N_smb_pos    !> SMB >= 0 North
52real                                 :: grad_S_smb_neg    !> SMB < 0 South
53real                                 :: grad_S_smb_pos    !> SMB >= 0 South
54real                                 :: coef_smb_unit     ! pour corriger l'unite
55
56real,dimension(nx,ny)                :: TA0               !> Temp annuelle sur S0
57
58
59! pour les lectures ncdf
60real*8, dimension(:,:,:),  pointer   :: tab3D             !< tableau 3d real pointer
61real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
62Character(len=10), dimension(3)      :: dimtestname       !> pour la sortie test netcdf
63integer                              :: ncidloc           !> pour la sortie test netcdf
64integer                              :: status            !> pour la sortie test netcdf
65
66integer                              :: massb_time        !< pour selectionner le type de calcul de smb
67! On a deux routines : celle avec un seul fichier (donnees observees) : massb_Ice2sea_fixe
68! Celle avec le bilan de masse sur plusieurs snapshots (annuels par ex.) entre lesquels on interpole.
69
70contains
71
72!--------------------------------------------------------------------------------
73!> SUBROUTINE: input_clim 
74!! Routine qui permet d'initialiser les variations temporelles des variables climatiques
75!>
76!--------------------------------------------------------------------------------
77
78subroutine input_clim 
79
80  character(len=100) :: smb_file           ! fichier smb
81  character(len=100) :: temp_annual_file   ! temperature annuelles
82  character(len=100) :: file_smb_snap      !> nom du fichier dans lequel sont les snapshots
83
84  integer            :: err                ! recuperation d'erreur
85  integer            :: i,j,k
86
87  namelist/clim_smb_T_gen/smb_file,coef_smb_unit,temp_annual_file
88
89428 format(A)
90  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
91  read(num_param,clim_smb_T_gen)
92
93  write(num_rep_42,428)'!________________________________________________________________' 
94  write(num_rep_42,428)'!  module climat_Grice2sea_years_mod lecture climat ref          '
95  write(num_rep_42,clim_smb_T_gen)
96  write(num_rep_42,428)'! smb_file          = fichier SMB (kg/m2/an)                     '
97  write(num_rep_42,428)'! coef_smb_unit     = coef passage m glace/an  (1/910 ou 1/918)  '
98  write(num_rep_42,428)'! temp_annual_file  = Temp moy annuelle  (°C)                    '
99  write(num_rep_42,428)'!________________________________________________________________'
100
101! smb : surface mass balance
102  smb_file  = trim(dirnameinp)//trim(smb_file)
103
104  call Read_Ncdf_var('z',smb_file,tab)
105
106!  call lect_input(3,'smb',1,bm,smb_file,trim(dirnameinp)//trim(runname)//'.nc')
107
108  bm(:,:)  = tab(:,:) * coef_smb_unit
109
110!  where ((H(:,:).lt.1.).and.(Bsoc(:,:).gt.0.))
111!     bm(:,:) = bm(:,:) - 5.                     ! pour faire un masque a l'exterieur du Groenland actuel
112!  end where
113
114!cdc test debug Hemin15 et Greeneem15
115!  where (bm(:,:).lt.-1000) bm(:,:)=0.
116
117
118  acc(:,:) = 0.
119  abl(:,:) = 0.
120
121  where (bm(:,:).gt.0.)
122     acc(:,:) = bm(:,:)   ! accumulation quand positif
123  elsewhere
124     abl(:,:) = - bm(:,:) ! ablation quand negatif
125  end where
126
127
128! surface temperature  Tann
129
130  temp_annual_file = trim(dirnameinp)//trim(temp_annual_file)
131
132  call Read_Ncdf_var('z',temp_annual_file,tab)
133  Tann(:,:)  = tab(:,:)
134!  call lect_input(3,'Tann',1,Tann,temp_annual_file,trim(dirnameinp)//trim(runname)//'.nc')
135
136!cdc test debug Hemin15
137!  where (Tann(:,:).lt.-100.) Tann(:,:)=10.
138
139
140  ta0(:,:)   = Tann(:,:)
141  Tjuly(:,:) = Tann(:,:)
142
143
144
145!lecture namelist
146
147  namelist/clim_snap/nb_snap,time_depart_snaps,ecart_snap,file_smb_snap,massb_time
148
149  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
150  read(num_param,clim_snap)
151
152! formats pour les ecritures dans 42
153  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
154  read(num_param,clim_snap)
155
156  write(num_rep_42,428)'!_______________________________________________________________________' 
157  write(num_rep_42,428)'!  module climat_Grice2sea_years_mod                                    '
158  write(num_rep_42,clim_snap)
159  write(num_rep_42,428)'! nb_snap           = nombre de snapshots                               '
160  write(num_rep_42,428)'! time_depart_snaps = debut du forçage                                  '
161  write(num_rep_42,428)'! ecart_snap        = ecart entre les snapshots                         '
162  write(num_rep_42,428)'! file_smb_snap     = fichier serie temp anomalie SMB de GCM            '
163  write(num_rep_42,428)'! massb_time        = 0:fixe, 1:interp snapshots, 2:snapsh+interp vert  '
164  write(num_rep_42,428)'!_______________________________________________________________________' 
165
166  if (massb_time == 1) then ! lecture des snapshots
167! allocation dynamique de time_snap
168     if (allocated(time_snap)) then
169        deallocate(time_snap,stat=err)
170        if (err/=0) then
171           print *,"Erreur à la desallocation de time_snap",err
172           stop
173        end if
174     end if
175
176     allocate(time_snap(nb_snap),stat=err)
177     if (err/=0) then
178        print *,"erreur a l'allocation du tableau time_snap ",err
179        print *,"nb_snap = ",nb_snap
180        stop
181     end if
182
183! remplissage de time_snap
184!write(6,*) 'time_snap'
185     do i=1,nb_snap
186        time_snap(i) = time_depart_snaps + ecart_snap * (i-1)
187!   write(6,*) i,time_snap(i)
188     end do
189
190! allocation dynamique de smb_snap
191     if (allocated(smb_snap)) then
192        deallocate(smb_snap,stat=err)
193        if (err/=0) then
194           print *,"Erreur à la desallocation de smb_snap",err
195           stop
196        end if
197     end if
198
199     allocate(smb_snap(nx,ny,nb_snap),stat=err)
200     if (err/=0) then
201        print *,"erreur a l'allocation du tableau smb_snap ",err
202        print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap
203        stop
204     end if
205
206! lecture de smb_snap
207     file_smb_snap  = trim(dirnameinp)//trim(file_smb_snap)
208     call Read_Ncdf_var('z',file_smb_snap,tab3D)
209     smb_snap (:,:,:) = Tab3D(:,:,:) * coef_smb_unit
210
211! ce sont des anomalies : ajoute les valeurs de reference
212     do k = 1,nb_snap
213        do j = 1,ny
214           do i = 1,nx
215              smb_snap (i,j,k) = smb_snap(i,j,k) + bm(i,j)
216           end do
217        end do
218     end do
219
220! copie la valeur de reference dans bm_0
221     bm_0(:,:) = bm(:,:)
222  endif
223
224! ecriture de verification
225!file_smb_snap = 'test_sortie_smb.nc'
226!file_smb_snap  = trim(dirnameout)//trim(file_smb_snap)
227
228!ncidloc = 501
229
230!call  lect_netcdf_type
231!write(6,*) 'ncdf_type'
232
233!if (ncdf_type.eq.32) then
234!   status  = nf90_create(TRIM(file_smb_snap),NF90_WRITE,ncidloc)    ! ouverture du fichier
235!else if (ncdf_type.eq.64) then
236!   status  = nf90_create(trim(file_smb_snap),and(nf90_write,nf90_64bit_offset),ncidloc) ! r2d2
237!end if
238
239!   status  = nf90_close(ncidloc)                                 ! fermeture
240
241! call write_ncdf_dim('x',trim(file_smb_snap),nx)               ! dimensions des tableaux
242! call write_ncdf_dim('y',trim(file_smb_snap),ny)
243! call write_ncdf_dim('time',trim(file_smb_snap),nb_snap)
244
245!Tab3D(:,:,:) = smb_snap (:,:,:)
246!dimtestname(1) = 'x'
247!dimtestname(2) = 'y'
248!dimtestname(3) = 'time'
249
250!call write_ncdf_var(trim('smb'),dimtestname,trim(file_smb_snap),tab3D,'double')
251   
252end subroutine input_clim
253
254!--------------------------------------------------------------------------------
255!> SUBROUTINE: init_forclim
256!! Routine qui permet d'initialiser les variables climatiques au cours du temps
257!>
258subroutine init_forclim
259
260  implicit none
261  namelist/lapse_rates/T_lapse_rate
262
263  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
264  read(num_param,lapse_rates)
265
266! formats pour les ecritures dans 42
267428 format(A)
268
269  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
270  read(num_param,lapse_rates)
271
272  write(num_rep_42,428)'!________________________________________________________________' 
273  write(num_rep_42,428)'!  module climat_Grice2sea_years_mod                             '
274  write(num_rep_42,lapse_rates)
275  write(num_rep_42,428)'!T_lapse_rate       = lapse rate temp annuelle                   '
276  write(num_rep_42,428)'!________________________________________________________________' 
277
278! appelle la routine de lecture des smb annuels
279  call input_clim
280  if (massb_time == 1) then ! lecture gradients smb
281     call init_grad_smb
282  endif
283
284  return
285end subroutine init_forclim
286
287!--------------------------------------------------------------------------------
288!> SUBROUTINE: forclim
289!! 
290!!  Routine qui permet le calcul climatique au cours du temps
291!!  @note Au temps considere (time) attribue les scalaires
292!!  @note   - tafor : forcage en temperature
293!!  @note   - sealevel : forcage niveau des mers
294!!  @note   - coefbmelt : forcage fusion basale ice shelves
295!>
296subroutine forclim               !  au temps considere (time)
297
298  implicit none
299
300  select case (massb_time)
301     case(0)
302! surface temperature et accumulation massb_Ice2sea_fixe
303! smb fixe et Tann avec lapse rate
304        Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) 
305        Ts(:,:)    = Tann(:,:)
306     case(1)
307        call massb_Ice2sea_RCM
308     case default
309        print *, 'Numero de massb invalide dans climat_Grice2sea_years_mod'
310        stop
311     end select
312end subroutine forclim
313
314
315subroutine massb_Ice2sea_RCM                ! calcule le mass balance
316
317  implicit none
318  integer             :: k_snap               ! pour calculer les indices de temps
319  real                :: time_bis             ! pour repliquer les dernieres annees
320  integer             :: k
321
322! calcul smb a partir fichier snapshots massb_Ice2sea_RCM
323! Calcule le mass balance d'apres un fichier de snapshots
324! avec la temperature parametree
325!     surface temperature et accumulation
326  Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) 
327  Ts(:,:)    = Tann(:,:)
328! calcule bm_time par interpolation entre deux snapshots
329! avant prend la valeur de reference
330! apres prend la derniere valeur
331! en general les snapshots vont de 2000 a 2200
332  if(time.lt.time_snap(1)) then              ! time avant le forcage
333     bm_time(:,:) = bm_0(:,:)
334     k_snap       = 1
335     S_ref(:,:)   = S(:,:)                   ! du coup sera la surface de reference avant le forcage
336     icum         = 0
337     i_moy        = 0
338     bm_ref_t(:,:)=  bm_0(:,:)               ! bilan de masse de reference
339     time_prec    = time
340  else if (time.ge.time_snap(nb_snap)) then  ! time apres le forcage
341     k_snap       = nb_snap
342     bm_time(:,:) =  smb_snap (:,:,k_snap)
343     if (abs(time-time_prec-1.).lt.dt) then   !
344        time_prec = time_prec + 1
345        icum      = 1
346     else
347        icum      = 0
348     endif
349  else                                       ! cas general
350     do k = 1 , nb_snap-1
351        if((time.ge.time_snap(k)).and.(time.lt.time_snap(k+1))) then ! entre k et k+1
352           bm_time(:,:) = smb_snap(:,:,k)+(smb_snap(:,:,k+1)-smb_snap(:,:,k)) *   &
353                (time-time_snap(k))/(time_snap(k+1)-time_snap(k))
354! exactement sur le snapshot et avec un ecart 1 an par rapport au precedent stockage
355!           write(6,*) 'time,tests',k,time,time-time_snap(k),time-time_prec-1.
356           if ((abs(time-time_snap(k)).le.dt).and.(abs(time-time_prec-1.).lt.dt)) then   
357              k_snap    = k
358              icum      = 1
359              time_prec = time_snap(k)     ! time_prec est le temps du precedent
360                                                 ! snapshot garde
361           else
362              icum    = 0
363           endif
364           exit
365        endif
366     end do
367  endif
368  call  grad_smb              !-----------------------------> A faire
369
370  if (massb_time == 1) then ! pas d'interpolation verticale
371     bm(:,:) = bm_time(:,:)
372  else if (massb_time == 2) then ! interpolation verticale
373     ! ajuste bm en fonction du temps et du gradient
374     bm(:,:) = bm_time(:,:) + grad_bm(:,:) *(S(:,:) - S_ref(:,:))
375     write(6,897) time, time_prec, icum, i_moy
376897  format('test temps smb   ',2(f0.3,1x),2(i0,1x))
377  endif
378
379! garde les 10 dernieres annees et calcule la moyenne
380  if (icum.eq.1) then                           ! stockage dans le tableau bm_ref_10
381     do k = 9,1,-1
382        bm_ref_10(:,:,k+1) = bm_ref_10(:,:,k)   ! on decale tous les elements
383     end do
384     bm_ref_10(:,:,k+1) = bm(:,:)               ! le plus recent est en position 1
385     i_moy              = i_moy +1              ! compte combien il y en a pour la moyenne
386     i_moy              = min(i_moy,10)
387     bm_ref_t(:,:)   = 0.
388     do k = 1,i_moy
389        bm_ref_t(:,:) = bm_ref_t(:,:) + bm_ref_10(:,:,k)
390     end do
391     bm_ref_t(:,:)    = bm_ref_t(:,:)/i_moy
392     write(6,898) time, time_prec, icum, i_moy
393898  format('cumul pour gradient  ',2(f0.3,1x),2(i0,1x))
394  end if
395end subroutine massb_Ice2sea_RCM
396 
397!------------------------------------------------------------------------------
398!> initialise le calcul du gradient vertical de smb
399subroutine init_grad_smb 
400
401  use module3D_phy
402  implicit none
403
404  character(len=120)       ::  file_grad_smb ! nom du fichier gradients de Tamsin
405  character(len=40)        ::  fin_ligne     ! fin de la ligne
406  real                     ::  grad
407
408  namelist/grad_smb/file_grad_smb
409
410! Dans lequel sont :
411! grad_N_smb_pos,grad_N_smb_neg,grad_S_smb_pos,grad_S_smb_neg,lim_lat
412
413  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
414  read(num_param,grad_smb)
415
416
417! formats pour les ecritures dans 42
418428 format(A)
419
420  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
421  read(num_param,grad_smb)
422
423  write(num_rep_42,428)'!________________________________________________________________' 
424  write(num_rep_42,428)'!  gradient smb climat_Grice2sea_years_mod                       '
425  write(num_rep_42,grad_smb)
426  write(num_rep_42,428)'!grad_smb       = fichier gradient SMB                           '
427  write(num_rep_42,428)'!________________________________________________________________' 
428
429
430  file_grad_smb=trim(dirnameinp)//'SMB-H-Feedback/gradients_18_07_2012/'//trim(file_grad_smb)
431  open(622,file=file_grad_smb)
432  do i=1,4
433     read(622,'(f9.4,A)') grad,fin_ligne
434     write(6,*) grad,fin_ligne
435
436     if (index(fin_ligne, "North").ne.0) then       ! North est dans la ligne
437        if (index(fin_ligne, "<").ne.0)  then       ! smb negatif
438           grad_N_smb_neg = grad
439        else if (index(fin_ligne, ">=").ne.0) then  ! smb positif
440           grad_N_smb_pos = grad
441        else
442           write(6,*) 'probleme lecture North fichier smb', file_grad_smb
443           STOP
444        end if
445
446     else if (index(fin_ligne, "South").ne.0) then   !South est dans la ligne
447        if (index(fin_ligne, "<").ne.0)  then        ! smb negatif
448           grad_S_smb_neg = grad
449        else if (index(fin_ligne, ">=").ne.0) then   ! smb positif
450           grad_S_smb_pos = grad
451        else
452           write(6,*) 'probleme lecture South fichier smb  ', file_grad_smb
453           STOP
454        end if
455
456     else
457        write(6,*) 'probleme lecture ni North ni South fichier smb  ', file_grad_smb
458        write(6,*) 'fin_ligne',fin_ligne,'   index North', index(fin_ligne, "North"),'   index South', index(fin_ligne, "South")
459        STOP
460     end if
461  end do
462
463  write(6,*) 'coefficients lus '
464  write(6,*) 'grad_N_smb_pos', grad_N_smb_pos
465  write(6,*) 'grad_N_smb_neg', grad_N_smb_neg
466  write(6,*) 'grad_S_smb_pos', grad_S_smb_pos
467  write(6,*) 'grad_S_smb_neg', grad_S_smb_neg
468
469  grad_N_smb_pos = coef_smb_unit * grad_N_smb_pos
470  grad_N_smb_neg = coef_smb_unit * grad_N_smb_neg
471  grad_S_smb_pos = coef_smb_unit * grad_S_smb_pos
472  grad_S_smb_neg = coef_smb_unit * grad_S_smb_neg
473
474  write(6,*) 'coefficients *coef_smb_unit '
475  write(6,*) 'grad_N_smb_pos', grad_N_smb_pos
476  write(6,*) 'grad_N_smb_neg', grad_N_smb_neg
477  write(6,*) 'grad_S_smb_pos', grad_S_smb_pos
478  write(6,*) 'grad_S_smb_neg', grad_S_smb_neg
479
480
481end subroutine init_grad_smb
482
483
484!------------------------------------------------------------------------------
485!> Calcule le gradient vertical de smb
486
487subroutine grad_smb 
488
489  use module3D_phy
490  implicit none
491
492  do j = 1,ny
493     do i = 1,nx
494        if (Ylat(i,j).gt.77.) then                     ! region nord
495           if (bm_ref_t(i,j).lt. 0.) then              ! smb negatif
496              grad_bm (i,j) = grad_N_smb_neg
497
498           else if (bm_ref_t(i,j).ge. 0.) then         ! smb positif
499              grad_bm (i,j) = grad_N_smb_pos 
500
501           end if
502
503        else if (Ylat(i,j).le.77.) then                ! region sud
504           if (bm_ref_t(i,j).lt. 0.) then              ! smb negatif
505              grad_bm (i,j) = grad_S_smb_neg
506
507           else if (bm_ref_t(i,j).ge. 0.) then         ! smb positif
508              grad_bm (i,j) = grad_S_smb_pos 
509
510           end if
511        end if
512     end do
513  end do
514
515end subroutine grad_smb
516
517end module  climat_Grice2sea_years_mod
Note: See TracBrowser for help on using the repository browser.