source: trunk/SOURCES/Ant16_files/climat_InitMIP_years_perturb_mod.f90 @ 288

Last change on this file since 288 was 288, checked in by dumas, 5 years ago

ISMIP forcing : smb and Tann without interpolation between years
File size: 15.2 KB
Line 
1!> \file climat_InitMIP_mod.f90
2! forcage avec BM
3!! Module ou les variations temporelles des variables climatiques
4!! sont directement imposees
5!<
6
7!> \namespace  climat_InitMIP_mod
8!! Module ou les variations temporelles des variables climatiques
9!! sont directement imposees en suivant le protocol initMIP
10!! \author afq
11!! \date 10 oct 17
12!! @note Used modules:
13!! @note   - use module3D_phy
14!<
15module climat_InitMIP_years_perturb_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,num_forc,dirforcage,dirnameinp,tafor,time,sealevel,sealevel_2d,coefbmshelf
19use netcdf
20use io_netcdf_grisli
21implicit none
22
23
24real :: T_lapse_rate              !< pour la temperature
25
26! anomalies: smb and Tann
27real,dimension(:,:,:),allocatable   :: smb_anom    !> bilan de masse des anomalies indices : nx,ny,tsnap
28real,dimension(:,:,:),allocatable   :: tann_anom   !> Tann anomalies indices : nx,ny,tsnap
29real,dimension(:),allocatable       :: time_snap   !> date des snapshots  indice : nb de nb_snap
30real,dimension(:,:,:),allocatable   :: grad_smb     !> gradient vertical de smb indices : nx,ny,tsnap
31real,dimension(:,:,:),allocatable   :: grad_tann   !> gradient vertical de smb indices : nx,ny,tsnap
32
33! declaration pour le bilan de masse et le bmelt
34real,dimension(nx,ny)    :: bm_0              !> bilan de masse initial
35
36real                     :: coef_smb_unit     ! pour corriger l'unite
37real                     :: coef_smb_anom_unit ! facteur correction unité anomalies de smb
38real,dimension(nx,ny)    :: TA0               !> Temp annuelle sur S0
39
40
41! aurel, pour climat type perturb:
42!integer nft             !> nombre de lignes a lire dans le fichier forcage climatique
43!real,dimension(:),allocatable :: tdate          !> time for climate forcing
44!real,dimension(:),allocatable :: tpert          !> temperature for climate forcing
45!real,dimension(:),allocatable :: spert          !> sea surface perturbation
46!real :: coefT                    !> pour modifier l'amplitude de la perturb. T
47!character(len=120) :: filforc    !> nom du fichier forcage
48!integer :: pertsmb               !> boolean: do we modify the smb?
49!real    :: rapsmb                !> if we modify the smb this is the equivalent of rappact
50
51
52! pour les lectures ncdf
53real*8, dimension(:), pointer        :: tab1D             !< tableau 1d real pointer
54real*8, dimension(:,:),   pointer    :: tab2D             !< tableau 2d real pointer
55real*8, dimension(:,:,:), pointer    :: tab3D             !< tableau 3d real pointer
56Character(len=10), dimension(3)      :: dimtestname       !> pour la sortie test netcdf
57integer                              :: ncidloc           !> pour la sortie test netcdf
58integer                              :: status            !> pour la sortie test netcdf
59
60integer                              :: massb_time        !< pour selectionner le type de calcul de smb
61
62integer                              :: nb_snap           !> nombre de snapshots
63
64real,dimension(nx,ny)                :: smb_anom_last20years  !> anomalie de smb des 20 dernieres annees du forcage
65real,dimension(nx,ny)                :: tann_anom_last20years !> anomalie de Tann des 20 dernieres annees du forcage
66
67contains
68
69!--------------------------------------------------------------------------------
70!> SUBROUTINE: input_clim 
71!! Routine qui permet d'initialiser les variations temporelles des variables climatiques
72!>
73!--------------------------------------------------------------------------------
74
75subroutine input_clim 
76
77  character(len=100) :: smb_file           ! fichier smb
78  character(len=100) :: temp_annual_file   ! temperature annuelles
79  character(len=100) :: file_smb_anom      !> nom du fichier dans lequel il y a l'anomalie smb
80  character(len=100) :: file_lapse         !> nom du fichier dans lequel il y a les lapserates (smb et tann)
81
82  integer            :: err                ! recuperation d'erreur
83  integer            :: i,k
84
85  !aurel for Tafor:
86  character(len=8)   :: control      !label to check clim. forc. file (filin) is usable
87  character(len=80)  :: filin
88 
89  real               :: time_depart_snaps  !> temps du debut premier snapshot
90
91  namelist/clim_smb_T_gen/smb_file,coef_smb_unit,temp_annual_file
92  namelist/smb_anom_initMIP/file_smb_anom,file_lapse,coef_smb_anom_unit,nb_snap,time_depart_snaps,massb_time
93 
94!!!!!  namelist/clim_snap/nb_snap,file_smb_snap,massb_time_snap
95
96428 format(A)
97  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
98  read(num_param,clim_smb_T_gen)
99
100  write(num_rep_42,428)'!________________________________________________________________' 
101  write(num_rep_42,428)'!  module climat_InitMIP_years_perturb_mod lecture climat ref          '
102  write(num_rep_42,clim_smb_T_gen)
103  write(num_rep_42,428)'! smb_file          = fichier SMB (kg/m2/an)                     '
104  write(num_rep_42,428)'! coef_smb_unit     = coef passage m glace/an  (1/910 ou 1/918)  '
105  write(num_rep_42,428)'! temp_annual_file  = Temp moy annuelle  (°C)                    '
106  write(num_rep_42,428)'!________________________________________________________________'
107
108 
109! smb : surface mass balance
110  smb_file  = trim(dirnameinp)//trim(smb_file)
111
112  call Read_Ncdf_var('smb',smb_file,tab2D)
113
114
115  bm(:,:)  = tab2D(:,:) * coef_smb_unit
116
117  !where ((H(:,:).lt.1.))
118  !   bm(:,:) = bm(:,:) - 15.                     ! pour faire un masque a l'exterieur du Groenland actuel
119  !end where
120
121  acc(:,:) = 0.
122  abl(:,:) = 0.
123
124  where (bm(:,:).gt.0.)
125     acc(:,:) = bm(:,:)   ! accumulation quand positif
126  elsewhere
127     abl(:,:) = - bm(:,:) ! ablation quand negatif
128  end where
129
130
131! surface temperature  Tann
132
133  temp_annual_file = trim(dirnameinp)//trim(temp_annual_file)
134
135  call Read_Ncdf_var('Tann',temp_annual_file,tab2D)
136  Tann(:,:)  = tab2D(:,:)
137
138  ta0(:,:)   = Tann(:,:)
139  Tjuly(:,:) = Tann(:,:)
140  bm_0(:,:) = bm(:,:)
141
142  sealevel=0.
143  tafor=0.
144  sealevel_2d(:,:) = sealevel
145
146  ! ______ Anomalies de SMB (simule asmb initMIP)
147  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
148  read(num_param,smb_anom_initMIP)
149
150  write(num_rep_42,428)'!____________________________________________________________________________' 
151  write(num_rep_42,428)'!  module climat_InitMIP_years_perturb_mod                                   '
152  write(num_rep_42,smb_anom_initMIP)
153  write(num_rep_42,428)'! file_smb_anom = fichier anomalie SMB de GCM                                '
154  write(num_rep_42,428)'! file_lapse    = fichier des lapse rate SMB et Tann                         '
155  write(num_rep_42,428)'! coef_smb_anom_unit = coef passage m glace/an  (1/910 ou 1/918)             '
156  write(num_rep_42,428)'! nb_snap       = nombre de snapshots                                        '
157  write(num_rep_42,428)'! time_depart_snaps = debut du forçage                                       '
158  write(num_rep_42,428)'! massb_time = 0:fixe, 1:anomalies, 2:interp snapshots, 3:snapsh+interp vert '
159  write(num_rep_42,428)'!____________________________________________________________________________'
160 
161  file_smb_anom  = trim(dirnameinp)//trim(file_smb_anom)
162  file_lapse     = trim(dirnameinp)//trim(file_lapse)
163
164  if (massb_time == 1) then ! fichier 2D de smb (exp initMIP asmb)
165! allocation dynamique de smb_snap
166      if (allocated(smb_anom)) then
167        deallocate(smb_anom,stat=err)
168        if (err/=0) then
169          print *,"Erreur à la desallocation de smb_anom",err
170          stop
171        end if
172      end if
173
174      allocate(smb_anom(nx,ny,1),stat=err)
175      if (err/=0) then
176        print *,"erreur a l'allocation du tableau smb_anom ",err
177        print *,"nx,ny,1 = ",nx,',',ny,',',1
178        stop
179      end if
180      call Read_Ncdf_var('asmb',file_smb_anom,tab2D)
181      smb_anom(:,:,1) = Tab2D(:,:) * coef_smb_anom_unit
182  elseif (massb_time >= 2) then ! fichier 3D de smb : lecture des snapshots
183! allocation dynamique de time_snap
184      if (allocated(time_snap)) then
185        deallocate(time_snap,stat=err)
186        if (err/=0) then
187          print *,"Erreur à la desallocation de time_snap",err
188          stop
189        end if
190      end if
191
192      allocate(time_snap(nb_snap),stat=err)
193      if (err/=0) then
194        print *,"erreur a l'allocation du tableau time_snap ",err
195        print *,"nb_snap = ",nb_snap
196        stop
197      end if
198! allocation dynamique de smb_snap
199      if (allocated(smb_anom)) then
200        deallocate(smb_anom,stat=err)
201        if (err/=0) then
202          print *,"Erreur à la desallocation de smb_anom",err
203          stop
204        end if
205      end if
206
207      allocate(smb_anom(nx,ny,nb_snap),stat=err)
208      if (err/=0) then
209        print *,"erreur a l'allocation du tableau smb_anom ",err
210        print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap
211        stop
212      end if
213
214! allocation dynamique de tann_snap
215      if (allocated(tann_anom)) then
216        deallocate(tann_anom,stat=err)
217        if (err/=0) then
218          print *,"Erreur à la desallocation de tann_anom",err
219          stop
220        end if
221      end if
222
223      allocate(tann_anom(nx,ny,nb_snap),stat=err)
224      if (err/=0) then
225        print *,"erreur a l'allocation du tableau tann_anom ",err
226        print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap
227        stop
228      end if
229 
230! lecture de smb_snap
231      call Read_Ncdf_var('smb_anomaly',file_smb_anom,tab3D)
232      smb_anom(:,:,:) = tab3D(:,:,:) * coef_smb_anom_unit
233 
234! lecture de tann_snap
235      call Read_Ncdf_var('ts_anomaly',file_smb_anom,tab3D)
236      tann_anom(:,:,:) = tab3D(:,:,:)
237
238! lecture de time_snap
239      call Read_Ncdf_var('time',file_smb_anom,tab1D)
240      time_snap(:) = tab1D(:)
241      time_snap(:) = time_snap(:) + time_depart_snaps
242   
243      if (massb_time == 3) then ! lecture lapse rate
244         
245         ! allocation dynamique de grad_smb
246         if (allocated(grad_smb)) then
247            deallocate(grad_smb,stat=err)
248            if (err/=0) then
249               print *,"Erreur à la desallocation de grad_smb",err
250               stop
251            end if
252         end if
253
254         allocate(grad_smb(nx,ny,nb_snap),stat=err)
255         if (err/=0) then
256            print *,"erreur a l'allocation du tableau grad_smb ",err
257            print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap
258            stop
259         end if
260
261     
262         ! allocation dynamique de smb_snap
263         if (allocated(grad_tann)) then
264            deallocate(grad_tann,stat=err)
265            if (err/=0) then
266               print *,"Erreur à la desallocation de grad_tann",err
267               stop
268            end if
269         end if
270
271         allocate(grad_tann(nx,ny,nb_snap),stat=err)
272         if (err/=0) then
273            print *,"erreur a l'allocation du tableau grad_tann ",err
274            print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap
275            stop
276         end if
277     
278         call Read_Ncdf_var('smb_gradient',file_lapse,tab3D)
279         grad_smb(:,:,:) = tab3D(:,:,:) * coef_smb_anom_unit
280         
281         call Read_Ncdf_var('ts_gradient',file_lapse,tab3D)
282         grad_tann(:,:,:) = tab3D(:,:,:)
283         
284      endif
285! calcul de la moyenne des 20 dernieres annees du forcage
286    do k=nb_snap,nb_snap-19,-1
287      smb_anom_last20years(:,:) = smb_anom_last20years(:,:) + smb_anom(:,:,k)
288      tann_anom_last20years(:,:) = tann_anom_last20years(:,:) + tann_anom(:,:,k)
289    enddo
290      smb_anom_last20years(:,:) = smb_anom_last20years(:,:)/20.
291      tann_anom_last20years(:,:) = tann_anom_last20years(:,:)/20.
292  endif     
293
294   
295end subroutine input_clim
296
297!--------------------------------------------------------------------------------
298!> SUBROUTINE: init_forclim
299!! Routine qui permet d'initialiser les variables climatiques au cours du temps
300!>
301subroutine init_forclim
302
303  implicit none
304 
305end subroutine init_forclim
306
307!--------------------------------------------------------------------------------
308!> SUBROUTINE: forclim
309!! 
310!!  Routine qui permet le calcul climatique au cours du temps
311!>
312
313subroutine forclim               !  au temps considere (time)
314
315  implicit none
316
317  integer i,ift
318  real             :: coefanomtime
319 
320 
321  if (massb_time == 0) then
322    bm(:,:) = bm_0(:,:)
323  elseif (massb_time == 1) then
324    coefanomtime = min ( real(time/40.) , 1. ) 
325    bm(:,:) = bm_0(:,:) + coefanomtime * smb_anom(:,:,1)
326  elseif (massb_time >= 2) then
327    call massb_ISMIP_RCM
328  end if
329 
330  Ts(:,:)    = Tann(:,:)
331
332       
333! coefmshelf est un coefficient qui fait varier bmgrz et bmshelf en fonction de tafor
334  coefbmshelf=1.
335     
336end subroutine forclim
337
338subroutine massb_ISMIP_RCM                ! calcule le mass balance
339
340  implicit none
341  integer             :: k                ! pour calculer les indices de temps
342  real,dimension(nx,ny) :: bm_time, tann_time ! pour calcul Bm et Tann
343  real,dimension(nx,ny) :: grad_smb_time, grad_tann_time ! pour calcul Bm et Tann avec correction verticale
344
345! calcul smb a partir fichier snapshots massb_ISMIP_RCM
346! Calcule le mass balance d'apres un fichier de snapshots
347
348! calcule bm_time par interpolation entre deux snapshots
349! avant prend la valeur de reference
350! apres prend la derniere valeur
351! en general les snapshots vont de 1995 a 2100 (time_snap de 0 à 105)
352  if(time.lt.time_snap(1)) then              ! time avant le forcage
353     bm_time(:,:) = bm_0(:,:)
354     tann_time(:,:) = ta0(:,:)
355  else if (time.ge.time_snap(nb_snap)) then  ! time apres le forcage : constant = derniere année nb_snap
356     bm_time(:,:) =  bm_0(:,:) + smb_anom (:,:,nb_snap)
357     tann_time(:,:) = ta0(:,:) + tann_anom(:,:,nb_snap)
358     if (massb_time == 3 ) then
359        grad_smb_time(:,:) =   grad_smb (:,:,1)
360        grad_tann_time(:,:) = grad_tann(:,:,1)
361     endif
362!~   else if (time.ge.time_snap(nb_snap)) then  ! time apres le forcage : constant moy 20 dernieres années
363!~      bm_time(:,:) =  bm_0(:,:) + smb_anom_last20years(:,:)
364!~      tann_time(:,:) = ta0(:,:) + tann_anom_last20years(:,:)
365  else                                       ! cas general
366     do k = 1 , nb_snap-1
367        if((time.ge.time_snap(k)).and.(time.lt.time_snap(k+1))) then ! entre k et k+1
368!cdc version avec interpolation lineaire entre 2 années       
369!~            bm_time(:,:) = bm_0(:,:) + smb_anom(:,:,k) + (smb_anom(:,:,k+1)-smb_anom(:,:,k)) *   &
370!~                 (time-time_snap(k))/(time_snap(k+1)-time_snap(k))
371!~            tann_time(:,:) = ta0(:,:) + tann_anom(:,:,k) + (tann_anom(:,:,k+1)-tann_anom(:,:,k)) *   &
372!~                 (time-time_snap(k))/(time_snap(k+1)-time_snap(k))
373!cdc version ISMIP : forcage annee n pendant un an sans interpolation.
374           bm_time(:,:) = bm_0(:,:) + smb_anom(:,:,k)
375           tann_time(:,:) = ta0(:,:) + tann_anom(:,:,k)
376           if (massb_time == 3 ) then
377              grad_smb_time(:,:) =   grad_smb(:,:,k) + (grad_smb(:,:,k+1)-grad_smb(:,:,k)) *   &
378                (time-time_snap(k))/(time_snap(k+1)-time_snap(k))
379              grad_tann_time(:,:) = grad_tann(:,:,k) + (grad_tann(:,:,k+1)-grad_tann(:,:,k)) *   &
380                (time-time_snap(k))/(time_snap(k+1)-time_snap(k))
381           endif
382           exit
383        endif
384     end do
385  endif
386
387  if (massb_time == 2) then ! pas d'interpolation verticale
388     bm(:,:) = bm_time(:,:)
389     Tann(:,:) = tann_time(:,:)
390  else if (massb_time == 3) then ! interpolation verticale
391     ! ajuste bm en fonction du temps et du gradient
392     bm(:,:) = bm_time(:,:) + grad_smb_time(:,:) * (S(:,:) - S0(:,:))
393     tann(:,:)= tann_time(:,:) + grad_tann_time(:,:) * (S(:,:) - S0(:,:))
394  endif
395
396end subroutine massb_ISMIP_RCM
397
398
399
400
401end module  climat_InitMIP_years_perturb_mod
Note: See TracBrowser for help on using the repository browser.