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