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 | !< |
---|
15 | module climat_Grice2sea_years_perturb_mod |
---|
16 | |
---|
17 | |
---|
18 | use 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,coefbmshelf |
---|
19 | !use lect_climref_Ice2sea |
---|
20 | use netcdf |
---|
21 | use io_netcdf_grisli |
---|
22 | implicit none |
---|
23 | |
---|
24 | |
---|
25 | real :: T_lapse_rate !< pour la temperature |
---|
26 | |
---|
27 | |
---|
28 | ! declarations specifiques year by year |
---|
29 | |
---|
30 | real,dimension(:),allocatable :: time_snap !> date des snapshots indice : nb de nb_snap) |
---|
31 | real,dimension(:,:,:),allocatable :: smb_snap !> bilan de masse des snapshots indices : nx,ny,nb_snap |
---|
32 | real :: ecart_snap !> ecart temporel entre les snapshots |
---|
33 | real :: time_depart_snaps !> temps du debut premier snapshot |
---|
34 | integer :: nb_snap !> nombre de snapshots |
---|
35 | |
---|
36 | ! declaration pour le bilan de masse |
---|
37 | real,dimension(nx,ny) :: bm_0 !> bilan de masse initial equ. S_0ref de Tamsin |
---|
38 | real,dimension(nx,ny) :: bm_time !> bilan de masse interpole entre 2 snapshots (~ S_t^RCM) |
---|
39 | |
---|
40 | ! calcul du gradient |
---|
41 | real,dimension(nx,ny) :: bm_ref_t !> bilan de masse de reference au temps t |
---|
42 | real,dimension(nx,ny,10) :: bm_ref_10 !> tableau des bm_ref pour moyenne 10 ans |
---|
43 | real :: time_prec !> temps du precedent snapshot |
---|
44 | integer :: icum !> pour le test stockage dans bm_ref_10 |
---|
45 | integer :: i_moy !> nombre de snapshots stockes |
---|
46 | |
---|
47 | real,dimension(nx,ny) :: grad_bm !> gradient du bilan de masse |
---|
48 | real,dimension(nx,ny) :: S_ref !> surface de reference |
---|
49 | |
---|
50 | real :: grad_N_smb_neg !> SMB < 0 North |
---|
51 | real :: grad_N_smb_pos !> SMB >= 0 North |
---|
52 | real :: grad_S_smb_neg !> SMB < 0 South |
---|
53 | real :: grad_S_smb_pos !> SMB >= 0 South |
---|
54 | real :: coef_smb_unit ! pour corriger l'unite |
---|
55 | |
---|
56 | real,dimension(nx,ny) :: TA0 !> Temp annuelle sur S0 |
---|
57 | |
---|
58 | |
---|
59 | ! aurel, pour climat type perturb: |
---|
60 | integer nft !> nombre de lignes a lire dans le fichier forcage climatique |
---|
61 | real,dimension(:),allocatable :: tdate !> time for climate forcing |
---|
62 | real,dimension(:),allocatable :: tpert !> temperature for climate forcing |
---|
63 | real,dimension(:),allocatable :: spert !> sea surface perturbation |
---|
64 | real,dimension(:),allocatable :: bpert !> basal melt index perturbation |
---|
65 | real :: coefT !> pour modifier l'amplitude de la perturb. T |
---|
66 | character(len=150) :: filforc !> nom du fichier forcage |
---|
67 | integer :: pertsmb !> boolean: do we modify the smb? |
---|
68 | real :: rapsmb !> if we modify the smb this is the equivalent of rappact |
---|
69 | integer :: pertbmb !> boolean: do we modify the bmelt under ice shelves? |
---|
70 | real :: coefbmb !> if we modify the bmelt, we read a perturbation (based on dD, insolation, etc.) and apply this coefficient |
---|
71 | |
---|
72 | |
---|
73 | |
---|
74 | ! pour les lectures ncdf |
---|
75 | real*8, dimension(:,:,:), pointer :: tab3D !< tableau 3d real pointer |
---|
76 | real*8, dimension(:,:), pointer :: tab !< tableau 2d real pointer |
---|
77 | Character(len=10), dimension(3) :: dimtestname !> pour la sortie test netcdf |
---|
78 | integer :: ncidloc !> pour la sortie test netcdf |
---|
79 | integer :: status !> pour la sortie test netcdf |
---|
80 | |
---|
81 | integer :: massb_time !< pour selectionner le type de calcul de smb |
---|
82 | ! On a deux routines : celle avec un seul fichier (donnees observees) : massb_Ice2sea_fixe |
---|
83 | ! Celle avec le bilan de masse sur plusieurs snapshots (annuels par ex.) entre lesquels on interpole. |
---|
84 | |
---|
85 | contains |
---|
86 | |
---|
87 | !-------------------------------------------------------------------------------- |
---|
88 | !> SUBROUTINE: input_clim |
---|
89 | !! Routine qui permet d'initialiser les variations temporelles des variables climatiques |
---|
90 | !> |
---|
91 | !-------------------------------------------------------------------------------- |
---|
92 | |
---|
93 | subroutine input_clim |
---|
94 | |
---|
95 | character(len=100) :: smb_file ! fichier smb |
---|
96 | character(len=100) :: temp_annual_file ! temperature annuelles |
---|
97 | character(len=100) :: file_smb_snap !> nom du fichier dans lequel sont les snapshots |
---|
98 | |
---|
99 | integer :: err ! recuperation d'erreur |
---|
100 | integer :: i,j,k |
---|
101 | |
---|
102 | !aurel for Tafor: |
---|
103 | character(len=8) :: control !label to check clim. forc. file (filin) is usable |
---|
104 | character(len=100):: filin |
---|
105 | |
---|
106 | namelist/clim_smb_T_gen/smb_file,coef_smb_unit,temp_annual_file |
---|
107 | namelist/clim_snap/nb_snap,time_depart_snaps,ecart_snap,file_smb_snap,massb_time |
---|
108 | |
---|
109 | rewind(num_param) ! pour revenir au debut du fichier param_list.dat |
---|
110 | read(num_param,clim_smb_T_gen) |
---|
111 | |
---|
112 | write(num_rep_42,'(A)')'! module climat_Grice2sea_years_perturb_mod ' |
---|
113 | write(num_rep_42,clim_smb_T_gen) |
---|
114 | write(num_rep_42,'(A)')'! smb_file = fichier SMB (kg/m2/an) ' |
---|
115 | write(num_rep_42,'(A)')'! coef_smb_unit = coef passage m glace/an (1/910 ou 1/918) ' |
---|
116 | write(num_rep_42,'(A)')'! temp_annual_file = Temp moy annuelle (°C) ' |
---|
117 | write(num_rep_42,'(A)')'!________________________________________________________________' |
---|
118 | |
---|
119 | |
---|
120 | ! smb : surface mass balance |
---|
121 | smb_file = trim(dirnameinp)//trim(smb_file) |
---|
122 | |
---|
123 | call Read_Ncdf_var('smb',smb_file,tab) |
---|
124 | |
---|
125 | ! call lect_input(3,'smb',1,bm,smb_file,trim(dirnameinp)//trim(runname)//'.nc') |
---|
126 | |
---|
127 | bm(:,:) = tab(:,:) * coef_smb_unit |
---|
128 | |
---|
129 | !where ((H(:,:).lt.1.).and.(Bsoc(:,:).gt.0.)) |
---|
130 | ! bm(:,:) = bm(:,:) - 5. ! pour faire un masque a l'exterieur du Groenland actuel |
---|
131 | !end where |
---|
132 | !where ((H(:,:).lt.1.).and.(Bsoc(:,:).gt.0.).and.(bm(:,:).gt.0.)) |
---|
133 | ! bm(:,:) = -0.1 |
---|
134 | !end where |
---|
135 | !cdc test debug Hemin15 et Greeneem15 |
---|
136 | ! where (bm(:,:).lt.-1000) bm(:,:)=0. |
---|
137 | where (bm(:,:).eq.0) bm(:,:)=-5. !afq |
---|
138 | |
---|
139 | acc(:,:) = 0. |
---|
140 | abl(:,:) = 0. |
---|
141 | |
---|
142 | where (bm(:,:).gt.0.) |
---|
143 | acc(:,:) = bm(:,:) ! accumulation quand positif |
---|
144 | elsewhere |
---|
145 | abl(:,:) = - bm(:,:) ! ablation quand negatif |
---|
146 | end where |
---|
147 | |
---|
148 | |
---|
149 | ! surface temperature Tann |
---|
150 | |
---|
151 | temp_annual_file = trim(dirnameinp)//trim(temp_annual_file) |
---|
152 | |
---|
153 | call Read_Ncdf_var('Tann',temp_annual_file,tab) |
---|
154 | Tann(:,:) = tab(:,:) |
---|
155 | ! call lect_input(3,'Tann',1,Tann,temp_annual_file,trim(dirnameinp)//trim(runname)//'.nc') |
---|
156 | |
---|
157 | !cdc test debug Hemin15 |
---|
158 | ! where (Tann(:,:).lt.-100.) Tann(:,:)=10. |
---|
159 | where (Tann(:,:).lt.-100.) Tann(:,:)=10. !afq |
---|
160 | |
---|
161 | ta0(:,:) = Tann(:,:) |
---|
162 | Tjuly(:,:) = Tann(:,:) |
---|
163 | |
---|
164 | rewind(num_param) ! pour revenir au debut du fichier param_list.dat |
---|
165 | read(num_param,clim_snap) |
---|
166 | |
---|
167 | write(num_rep_42,'(A)')'! module climat_Grice2sea_years_perturb_mod ' |
---|
168 | write(num_rep_42,clim_snap) |
---|
169 | write(num_rep_42,'(A)')'! nb_snap = nombre de snapshots ' |
---|
170 | write(num_rep_42,'(A)')'! time_depart_snaps = debut du forçage ' |
---|
171 | write(num_rep_42,'(A)')'! ecart_snap = ecart entre les snapshots ' |
---|
172 | write(num_rep_42,'(A)')'! file_smb_snap = fichier serie temp anomalie SMB de GCM ' |
---|
173 | write(num_rep_42,'(A)')'! massb_time = 0:fixe, 1:interp snapshots, 2:snapsh+interp vert ' |
---|
174 | write(num_rep_42,'(A)')'!_______________________________________________________________________' |
---|
175 | |
---|
176 | if (massb_time == 1) then ! lecture des snapshots |
---|
177 | ! allocation dynamique de time_snap |
---|
178 | if (allocated(time_snap)) then |
---|
179 | deallocate(time_snap,stat=err) |
---|
180 | if (err/=0) then |
---|
181 | print *,"Erreur à la desallocation de time_snap",err |
---|
182 | stop |
---|
183 | end if |
---|
184 | end if |
---|
185 | |
---|
186 | allocate(time_snap(nb_snap),stat=err) |
---|
187 | if (err/=0) then |
---|
188 | print *,"erreur a l'allocation du tableau time_snap ",err |
---|
189 | print *,"nb_snap = ",nb_snap |
---|
190 | stop |
---|
191 | end if |
---|
192 | |
---|
193 | ! remplissage de time_snap |
---|
194 | !write(6,*) 'time_snap' |
---|
195 | do i=1,nb_snap |
---|
196 | time_snap(i) = time_depart_snaps + ecart_snap * (i-1) |
---|
197 | ! write(6,*) i,time_snap(i) |
---|
198 | end do |
---|
199 | |
---|
200 | ! allocation dynamique de smb_snap |
---|
201 | if (allocated(smb_snap)) then |
---|
202 | deallocate(smb_snap,stat=err) |
---|
203 | if (err/=0) then |
---|
204 | print *,"Erreur à la desallocation de smb_snap",err |
---|
205 | stop |
---|
206 | end if |
---|
207 | end if |
---|
208 | |
---|
209 | allocate(smb_snap(nx,ny,nb_snap),stat=err) |
---|
210 | if (err/=0) then |
---|
211 | print *,"erreur a l'allocation du tableau smb_snap ",err |
---|
212 | print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap |
---|
213 | stop |
---|
214 | end if |
---|
215 | |
---|
216 | ! lecture de smb_snap |
---|
217 | file_smb_snap = trim(dirnameinp)//trim(file_smb_snap) |
---|
218 | call Read_Ncdf_var('z',file_smb_snap,tab3D) |
---|
219 | smb_snap (:,:,:) = Tab3D(:,:,:) * coef_smb_unit |
---|
220 | |
---|
221 | ! ce sont des anomalies : ajoute les valeurs de reference |
---|
222 | do k = 1,nb_snap |
---|
223 | do j = 1,ny |
---|
224 | do i = 1,nx |
---|
225 | smb_snap (i,j,k) = smb_snap(i,j,k) + bm(i,j) |
---|
226 | end do |
---|
227 | end do |
---|
228 | end do |
---|
229 | |
---|
230 | ! copie la valeur de reference dans bm_0 |
---|
231 | ! bm_0(:,:) = bm(:,:) |
---|
232 | endif |
---|
233 | ! copie la valeur de reference dans bm_0 |
---|
234 | bm_0(:,:) = bm(:,:) !afq marion dufresne: we do it even if we don't read any snapshot |
---|
235 | |
---|
236 | filin=trim(dirforcage)//trim(filforc) |
---|
237 | |
---|
238 | open(num_forc,file=filin,status='old') |
---|
239 | |
---|
240 | read(num_forc,*) control,nft |
---|
241 | |
---|
242 | ! Determination of file size (line nb), allocation of perturbation array |
---|
243 | |
---|
244 | if (control.ne.'nb_lines') then |
---|
245 | write(6,*) filin,'indiquer le nb de ligne en debut de fichier:' |
---|
246 | write(6,*) 'le nb de lignes et le label de control nb_lines' |
---|
247 | stop |
---|
248 | endif |
---|
249 | |
---|
250 | ! Dimensionnement des tableaux tdate, .... |
---|
251 | |
---|
252 | if (.not.allocated(tdate)) then |
---|
253 | allocate(tdate(nft),stat=err) |
---|
254 | if (err/=0) then |
---|
255 | print *,"erreur a l'allocation du tableau Tdate",err |
---|
256 | stop 4 |
---|
257 | end if |
---|
258 | end if |
---|
259 | |
---|
260 | if (.not.allocated(spert)) then |
---|
261 | allocate(spert(nft),stat=err) |
---|
262 | if (err/=0) then |
---|
263 | print *,"erreur a l'allocation du tableau Spert",err |
---|
264 | stop 4 |
---|
265 | end if |
---|
266 | end if |
---|
267 | |
---|
268 | if (.not.allocated(tpert)) then |
---|
269 | allocate(tpert(nft),stat=err) |
---|
270 | if (err/=0) then |
---|
271 | print *,"erreur a l'allocation du tableau Tpert",err |
---|
272 | stop 4 |
---|
273 | end if |
---|
274 | end if |
---|
275 | |
---|
276 | if (.not.allocated(bpert)) then |
---|
277 | allocate(bpert(nft),stat=err) |
---|
278 | if (err/=0) then |
---|
279 | print *,"erreur a l'allocation du tableau Bpert",err |
---|
280 | stop 4 |
---|
281 | end if |
---|
282 | end if |
---|
283 | |
---|
284 | do i=1,nft |
---|
285 | if (pertbmb==1) then |
---|
286 | read(num_forc,*) tdate(i),spert(i),tpert(i),bpert(i) |
---|
287 | else |
---|
288 | read(num_forc,*) tdate(i),spert(i),tpert(i) |
---|
289 | bpert(i)=0. |
---|
290 | endif |
---|
291 | end do |
---|
292 | close(num_forc) |
---|
293 | |
---|
294 | tpert(:)=tpert(:)*coefT |
---|
295 | bpert(:)=max(1.+bpert(:)*coefbmb,0.01) !freezing is not allowed |
---|
296 | |
---|
297 | if (massb_time == 1) then ! lecture gradients smb |
---|
298 | call init_grad_smb |
---|
299 | endif |
---|
300 | |
---|
301 | end subroutine input_clim |
---|
302 | |
---|
303 | !-------------------------------------------------------------------------------- |
---|
304 | !> SUBROUTINE: init_forclim |
---|
305 | !! Routine qui permet d'initialiser les variables climatiques au cours du temps |
---|
306 | !> |
---|
307 | subroutine init_forclim |
---|
308 | |
---|
309 | implicit none |
---|
310 | namelist/lapse_rates/T_lapse_rate |
---|
311 | namelist/clim_pert_massb/coefT,filforc,pertsmb,rapsmb,pertbmb,coefbmb |
---|
312 | |
---|
313 | rewind(num_param) ! pour revenir au debut du fichier param_list.dat |
---|
314 | read(num_param,lapse_rates) |
---|
315 | |
---|
316 | write(num_rep_42,'(A)')'! module climat_Grice2sea_years_perturb_mod ' |
---|
317 | write(num_rep_42,lapse_rates) |
---|
318 | write(num_rep_42,'(A)')'!T_lapse_rate = lapse rate temp annuelle ' |
---|
319 | write(num_rep_42,'(A)')'!________________________________________________________________' |
---|
320 | |
---|
321 | |
---|
322 | rewind(num_param) |
---|
323 | read(num_param,clim_pert_massb) |
---|
324 | |
---|
325 | write(num_rep_42,'(A)') '! module climat_Grice2sea_years_perturb_mod ' |
---|
326 | write(num_rep_42,clim_pert_massb) |
---|
327 | write(num_rep_42,'(A)') '! coefT : coef amplitude perturbation T' |
---|
328 | write(num_rep_42,'(A)') '! filforc : fichier de forcage temporel' |
---|
329 | write(num_rep_42,'(A)') '! pertsmb : 0 SMB fixe, 1 SMB exp temperature' |
---|
330 | write(num_rep_42,'(A)') '! rapsmb : coef variation SMB' |
---|
331 | write(num_rep_42,'(A)') '! pertbmb : 0 bmelt fixe, 1 bmelt variable' |
---|
332 | write(num_rep_42,'(A)') '! coefbmb : coef variation bmelt' |
---|
333 | write(num_rep_42,'(A)') '!________________________________________________________________' |
---|
334 | |
---|
335 | return |
---|
336 | end subroutine init_forclim |
---|
337 | |
---|
338 | !-------------------------------------------------------------------------------- |
---|
339 | !> SUBROUTINE: forclim |
---|
340 | !! |
---|
341 | !! Routine qui permet le calcul climatique au cours du temps |
---|
342 | !! @note Au temps considere (time) attribue les scalaires |
---|
343 | !! @note - tafor : forcage en temperature |
---|
344 | !! @note - sealevel : forcage niveau des mers |
---|
345 | !! @note - coefbmelt : forcage fusion basale ice shelves |
---|
346 | !> |
---|
347 | subroutine forclim ! au temps considere (time) |
---|
348 | |
---|
349 | implicit none |
---|
350 | |
---|
351 | integer i,ift |
---|
352 | |
---|
353 | select case (massb_time) |
---|
354 | case(0) |
---|
355 | ! smb fixe et Tann avec lapse rate + Tafor |
---|
356 | |
---|
357 | if(time.lt.tdate(1)) then |
---|
358 | tafor=tpert(1) |
---|
359 | sealevel=spert(1) |
---|
360 | coefbmshelf=bpert(1) |
---|
361 | ift=1 |
---|
362 | |
---|
363 | else if (time.ge.tdate(nft)) then |
---|
364 | tafor=tpert(nft) |
---|
365 | sealevel=spert(nft) |
---|
366 | coefbmshelf=bpert(nft) |
---|
367 | ift=nft |
---|
368 | |
---|
369 | else |
---|
370 | do i=1,nft-1 |
---|
371 | if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then ! entre i et i+1 : cas general |
---|
372 | tafor=tpert(i)+(tpert(i+1)-tpert(i))* & |
---|
373 | (time-tdate(i))/(tdate(i+1)-tdate(i)) |
---|
374 | sealevel=spert(i)+(spert(i+1)-spert(i))* & |
---|
375 | (time-tdate(i))/(tdate(i+1)-tdate(i)) |
---|
376 | coefbmshelf=bpert(i)+(bpert(i+1)-bpert(i))* & |
---|
377 | (time-tdate(i))/(tdate(i+1)-tdate(i)) |
---|
378 | ift=i |
---|
379 | goto 100 |
---|
380 | endif |
---|
381 | end do |
---|
382 | endif |
---|
383 | 100 continue |
---|
384 | |
---|
385 | Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor |
---|
386 | Ts(:,:) = Tann(:,:) |
---|
387 | |
---|
388 | ! aurel marion dufresne: we might want to decrease the SMB during glacials..? |
---|
389 | if (pertsmb.eq.1) then |
---|
390 | bm(:,:) = bm_0(:,:) * exp( rapsmb *(Tann(:,:)-Ta0(:,:))) |
---|
391 | if (Tafor.lt.0.) then |
---|
392 | where(bm(:,:).lt.0.) bm(:,:)=min(bm(:,:)-Tafor*0.05,1.) !10 degrees less give 0.5 meter more ? |
---|
393 | end if |
---|
394 | end if |
---|
395 | |
---|
396 | !afq -- as in Pollard and Deconto, melt is more efficient at high temperature... -> |
---|
397 | if (coefbmshelf.gt.1.) then |
---|
398 | coefbmshelf = 1. + ( coefbmshelf - 1. ) * 1.5 !afq.... tuning for the degla?... |
---|
399 | end if |
---|
400 | |
---|
401 | !afq --old: ! coefmshelf est un coefficient qui fait varier bmgrz et bmshelf en fonction de tafor |
---|
402 | !afq --old: coefbmshelf=(1.+tafor/10.) ! coefbmshelf=0 pour tafor=-10deg |
---|
403 | !afq --old: coefbmshelf=max(coefbmshelf,0.) |
---|
404 | !afq --old: coefbmshelf=min(coefbmshelf,2.) |
---|
405 | |
---|
406 | ! coefbmshelf=1. !afq for calibration Greenland |
---|
407 | |
---|
408 | case(1) |
---|
409 | call massb_Ice2sea_RCM |
---|
410 | case default |
---|
411 | print *, 'Numero de massb invalide dans climat_Grice2sea_years_mod' |
---|
412 | stop |
---|
413 | end select |
---|
414 | end subroutine forclim |
---|
415 | |
---|
416 | |
---|
417 | subroutine massb_Ice2sea_RCM ! calcule le mass balance |
---|
418 | |
---|
419 | implicit none |
---|
420 | integer :: k_snap ! pour calculer les indices de temps |
---|
421 | integer :: k |
---|
422 | |
---|
423 | ! calcul smb a partir fichier snapshots massb_Ice2sea_RCM |
---|
424 | ! Calcule le mass balance d'apres un fichier de snapshots |
---|
425 | ! avec la temperature parametree |
---|
426 | ! surface temperature et accumulation |
---|
427 | Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) |
---|
428 | Ts(:,:) = Tann(:,:) |
---|
429 | ! calcule bm_time par interpolation entre deux snapshots |
---|
430 | ! avant prend la valeur de reference |
---|
431 | ! apres prend la derniere valeur |
---|
432 | ! en general les snapshots vont de 2000 a 2200 |
---|
433 | if(time.lt.time_snap(1)) then ! time avant le forcage |
---|
434 | bm_time(:,:) = bm_0(:,:) |
---|
435 | k_snap = 1 |
---|
436 | S_ref(:,:) = S(:,:) ! du coup sera la surface de reference avant le forcage |
---|
437 | icum = 0 |
---|
438 | i_moy = 0 |
---|
439 | bm_ref_t(:,:)= bm_0(:,:) ! bilan de masse de reference |
---|
440 | time_prec = time |
---|
441 | else if (time.ge.time_snap(nb_snap)) then ! time apres le forcage |
---|
442 | k_snap = nb_snap |
---|
443 | bm_time(:,:) = smb_snap (:,:,k_snap) |
---|
444 | if (abs(time-time_prec-1.).lt.dt) then ! |
---|
445 | time_prec = time_prec + 1 |
---|
446 | icum = 1 |
---|
447 | else |
---|
448 | icum = 0 |
---|
449 | endif |
---|
450 | else ! cas general |
---|
451 | do k = 1 , nb_snap-1 |
---|
452 | if((time.ge.time_snap(k)).and.(time.lt.time_snap(k+1))) then ! entre k et k+1 |
---|
453 | bm_time(:,:) = smb_snap(:,:,k)+(smb_snap(:,:,k+1)-smb_snap(:,:,k)) * & |
---|
454 | (time-time_snap(k))/(time_snap(k+1)-time_snap(k)) |
---|
455 | ! exactement sur le snapshot et avec un ecart 1 an par rapport au precedent stockage |
---|
456 | ! write(6,*) 'time,tests',k,time,time-time_snap(k),time-time_prec-1. |
---|
457 | if ((abs(time-time_snap(k)).le.dt).and.(abs(time-time_prec-1.).lt.dt)) then |
---|
458 | k_snap = k |
---|
459 | icum = 1 |
---|
460 | time_prec = time_snap(k) ! time_prec est le temps du precedent |
---|
461 | ! snapshot garde |
---|
462 | else |
---|
463 | icum = 0 |
---|
464 | endif |
---|
465 | exit |
---|
466 | endif |
---|
467 | end do |
---|
468 | endif |
---|
469 | call grad_smb !-----------------------------> A faire |
---|
470 | |
---|
471 | if (massb_time == 1) then ! pas d'interpolation verticale |
---|
472 | bm(:,:) = bm_time(:,:) |
---|
473 | else if (massb_time == 2) then ! interpolation verticale |
---|
474 | ! ajuste bm en fonction du temps et du gradient |
---|
475 | bm(:,:) = bm_time(:,:) + grad_bm(:,:) *(S(:,:) - S_ref(:,:)) |
---|
476 | write(6,897) time, time_prec, icum, i_moy |
---|
477 | 897 format('test temps smb ',2(f0.3,1x),2(i0,1x)) |
---|
478 | endif |
---|
479 | |
---|
480 | ! garde les 10 dernieres annees et calcule la moyenne |
---|
481 | if (icum.eq.1) then ! stockage dans le tableau bm_ref_10 |
---|
482 | do k = 9,1,-1 |
---|
483 | bm_ref_10(:,:,k+1) = bm_ref_10(:,:,k) ! on decale tous les elements |
---|
484 | end do |
---|
485 | bm_ref_10(:,:,k+1) = bm(:,:) ! le plus recent est en position 1 |
---|
486 | i_moy = i_moy +1 ! compte combien il y en a pour la moyenne |
---|
487 | i_moy = min(i_moy,10) |
---|
488 | bm_ref_t(:,:) = 0. |
---|
489 | do k = 1,i_moy |
---|
490 | bm_ref_t(:,:) = bm_ref_t(:,:) + bm_ref_10(:,:,k) |
---|
491 | end do |
---|
492 | bm_ref_t(:,:) = bm_ref_t(:,:)/i_moy |
---|
493 | write(6,898) time, time_prec, icum, i_moy |
---|
494 | 898 format('cumul pour gradient ',2(f0.3,1x),2(i0,1x)) |
---|
495 | end if |
---|
496 | end subroutine massb_Ice2sea_RCM |
---|
497 | |
---|
498 | !------------------------------------------------------------------------------ |
---|
499 | !> initialise le calcul du gradient vertical de smb |
---|
500 | subroutine init_grad_smb |
---|
501 | |
---|
502 | use module3D_phy |
---|
503 | implicit none |
---|
504 | |
---|
505 | character(len=120) :: file_grad_smb ! nom du fichier gradients de Tamsin |
---|
506 | character(len=40) :: fin_ligne ! fin de la ligne |
---|
507 | real :: grad |
---|
508 | |
---|
509 | namelist/grad_smb/file_grad_smb |
---|
510 | |
---|
511 | ! Dans lequel sont : |
---|
512 | ! grad_N_smb_pos,grad_N_smb_neg,grad_S_smb_pos,grad_S_smb_neg,lim_lat |
---|
513 | |
---|
514 | rewind(num_param) ! pour revenir au debut du fichier param_list.dat |
---|
515 | read(num_param,grad_smb) |
---|
516 | |
---|
517 | |
---|
518 | ! formats pour les ecritures dans 42 |
---|
519 | 428 format(A) |
---|
520 | |
---|
521 | rewind(num_param) ! pour revenir au debut du fichier param_list.dat |
---|
522 | read(num_param,grad_smb) |
---|
523 | |
---|
524 | write(num_rep_42,428)'!________________________________________________________________' |
---|
525 | write(num_rep_42,428)'! gradient smb climat_Grice2sea_years_mod ' |
---|
526 | write(num_rep_42,grad_smb) |
---|
527 | write(num_rep_42,428)'!grad_smb = fichier gradient SMB ' |
---|
528 | write(num_rep_42,428)'!________________________________________________________________' |
---|
529 | |
---|
530 | |
---|
531 | file_grad_smb=trim(dirnameinp)//'SMB-H-Feedback/gradients_18_07_2012/'//trim(file_grad_smb) |
---|
532 | open(622,file=file_grad_smb) |
---|
533 | do i=1,4 |
---|
534 | read(622,'(f9.4,A)') grad,fin_ligne |
---|
535 | write(6,*) grad,fin_ligne |
---|
536 | |
---|
537 | if (index(fin_ligne, "North").ne.0) then ! North est dans la ligne |
---|
538 | if (index(fin_ligne, "<").ne.0) then ! smb negatif |
---|
539 | grad_N_smb_neg = grad |
---|
540 | else if (index(fin_ligne, ">=").ne.0) then ! smb positif |
---|
541 | grad_N_smb_pos = grad |
---|
542 | else |
---|
543 | write(6,*) 'probleme lecture North fichier smb', file_grad_smb |
---|
544 | STOP |
---|
545 | end if |
---|
546 | |
---|
547 | else if (index(fin_ligne, "South").ne.0) then !South est dans la ligne |
---|
548 | if (index(fin_ligne, "<").ne.0) then ! smb negatif |
---|
549 | grad_S_smb_neg = grad |
---|
550 | else if (index(fin_ligne, ">=").ne.0) then ! smb positif |
---|
551 | grad_S_smb_pos = grad |
---|
552 | else |
---|
553 | write(6,*) 'probleme lecture South fichier smb ', file_grad_smb |
---|
554 | STOP |
---|
555 | end if |
---|
556 | |
---|
557 | else |
---|
558 | write(6,*) 'probleme lecture ni North ni South fichier smb ', file_grad_smb |
---|
559 | write(6,*) 'fin_ligne',fin_ligne,' index North', index(fin_ligne, "North"),' index South', index(fin_ligne, "South") |
---|
560 | STOP |
---|
561 | end if |
---|
562 | end do |
---|
563 | |
---|
564 | write(6,*) 'coefficients lus ' |
---|
565 | write(6,*) 'grad_N_smb_pos', grad_N_smb_pos |
---|
566 | write(6,*) 'grad_N_smb_neg', grad_N_smb_neg |
---|
567 | write(6,*) 'grad_S_smb_pos', grad_S_smb_pos |
---|
568 | write(6,*) 'grad_S_smb_neg', grad_S_smb_neg |
---|
569 | |
---|
570 | grad_N_smb_pos = coef_smb_unit * grad_N_smb_pos |
---|
571 | grad_N_smb_neg = coef_smb_unit * grad_N_smb_neg |
---|
572 | grad_S_smb_pos = coef_smb_unit * grad_S_smb_pos |
---|
573 | grad_S_smb_neg = coef_smb_unit * grad_S_smb_neg |
---|
574 | |
---|
575 | write(6,*) 'coefficients *coef_smb_unit ' |
---|
576 | write(6,*) 'grad_N_smb_pos', grad_N_smb_pos |
---|
577 | write(6,*) 'grad_N_smb_neg', grad_N_smb_neg |
---|
578 | write(6,*) 'grad_S_smb_pos', grad_S_smb_pos |
---|
579 | write(6,*) 'grad_S_smb_neg', grad_S_smb_neg |
---|
580 | |
---|
581 | |
---|
582 | end subroutine init_grad_smb |
---|
583 | |
---|
584 | |
---|
585 | !------------------------------------------------------------------------------ |
---|
586 | !> Calcule le gradient vertical de smb |
---|
587 | |
---|
588 | subroutine grad_smb |
---|
589 | |
---|
590 | use module3D_phy |
---|
591 | implicit none |
---|
592 | |
---|
593 | do j = 1,ny |
---|
594 | do i = 1,nx |
---|
595 | if (Ylat(i,j).gt.77.) then ! region nord |
---|
596 | if (bm_ref_t(i,j).lt. 0.) then ! smb negatif |
---|
597 | grad_bm (i,j) = grad_N_smb_neg |
---|
598 | |
---|
599 | else if (bm_ref_t(i,j).ge. 0.) then ! smb positif |
---|
600 | grad_bm (i,j) = grad_N_smb_pos |
---|
601 | |
---|
602 | end if |
---|
603 | |
---|
604 | else if (Ylat(i,j).le.77.) then ! region sud |
---|
605 | if (bm_ref_t(i,j).lt. 0.) then ! smb negatif |
---|
606 | grad_bm (i,j) = grad_S_smb_neg |
---|
607 | |
---|
608 | else if (bm_ref_t(i,j).ge. 0.) then ! smb positif |
---|
609 | grad_bm (i,j) = grad_S_smb_pos |
---|
610 | |
---|
611 | end if |
---|
612 | end if |
---|
613 | end do |
---|
614 | end do |
---|
615 | |
---|
616 | end subroutine grad_smb |
---|
617 | |
---|
618 | end module climat_Grice2sea_years_perturb_mod |
---|