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

Last change on this file since 33 was 33, checked in by roche, 8 years ago

Small modification of the code to allow compilation in gfortran (moved a namelist declaration)

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