source: trunk/SOURCES/climat-forcage-insolation_mod.f90 @ 334

Last change on this file since 334 was 237, checked in by aquiquet, 5 years ago

Sealevel is now treated as a 2D variable (sealevel_2d while sealevel remains the eustatic sea level), results should remain identical as sealevel_2d is equal to sealevel in this revision.
File size: 43.9 KB
Line 
1!DERNIERE VERSION DU SCRIPT D'INTERPOLATION (May 2014) v2.1
2!Interpolation NL/lineaire rayon 200/400 km + orbital forcing + CO2 decrease interpo lin ou log
3
4module climat_forcage_insolation_mod                   
5
6! forcage avec champs mensuels
7! ce qui a trait aux lacs proglaciaires passe dans un module separe.
8! les tests sur geoplace sont enleves et sont remplacés par les
9! lectures de fichers appropries.
10! fonctionne avec un index en co2 et des snapshots a differents tx de co2
11! nouvelle version avec liste des variables utilisées par le module
12! C. Dumas 07/2015
13
14use module3D_phy,only:nx,ny,S,H,flot,Tann,Tjuly,Tmois,acc,coefbmshelf,sealevel_2d,ro,num_param,num_rep_42,dirnameinp,time,pi,dx,xlong,ylat,bsoc         
15use netcdf
16use io_netcdf_grisli
17       
18implicit none
19
20! 1=decalaration variables
21!-------------------------
22
23integer :: nft          !   NFT est le nombre de lignes a lire dans le fichier contenant le forcage climatique
24integer :: np           ! nbr de points englaces poses
25                       
26integer,parameter :: mois=12
27integer,parameter :: ntr=2      ! nb de snapshots selon les paramètres orbitaux
28integer,parameter :: gtr=4      ! nb de snapshots selon l'état de la calotte
29integer,parameter :: ctr=4      ! nb de snapshots selon le CO2
30real,dimension(ctr) :: Seuil_haut
31real,dimension(ctr) :: Seuil_bas
32integer :: CO2SnapMin,CO2SnapMax
33character(len=12) :: orbital_interpolation='La2004'   !Choose between La2004, constant and sinusoid
34character(len=12) :: ice_int='yes'         !Choose between yes and no
35character(len=12) :: interpolation_ice='nonlinear'   !Choose between linear and nonlinear
36integer :: rayon_interpolation_ice=60  !(In km) Choose a multiple of nx.
37
38!real :: poids_nb_points=0.5
39!real :: poids_proximity=0.5
40character(len=15) :: CO2_function='variable'   !Choose between variable and constant
41integer :: nb_slope_dis_CO2  ! nbr de ligne fichier CO2
42real :: CO2_value=980.
43character(len=12) :: interpolation_CO2='logarithmic'   !Choose between linear and logarithmic
44
45!integer :: itr        ! index of snapshot                     
46real :: coefsin                                          ! coef pour param orbitaux
47double precision :: factsin
48real :: orbite
49real :: coefCO2,ccoefCO2
50real :: factCO2
51real :: surf_test
52real :: retroicecoeff
53real,dimension(ntr) :: summorb  !Contient la moyenne DJF ou JJA de l'insolation pour
54! chacun des paramètres orbitaux. On a WSO = 485 W.m-² et CSO = 382.2 W.m-²
55real,dimension(ctr,gtr) :: palier_ice
56real,dimension(ctr) :: palier_CO2
57
58real,dimension(nx,ny,mois,ntr) :: Tm                    ! temperature mensuelle de chaque tranche
59real,dimension(nx,ny,mois,ntr) :: Pm                    ! precipitation mensuelle
60
61real,dimension(nx,ny,mois,ctr,gtr,ntr) :: Tm_fin        ! tableau d'interpolation stylé
62real,dimension(nx,ny,mois,ctr,gtr,ntr) :: Pm_fin        ! tableau d'interpolation stylé
63
64!interpolation sur param orbitaux
65
66!Définition des variables propres au calcul de la fonction d'Insolation
67!Nom du fichier contenant les points de temps et d'insolation récupérés du site
68!de Laskar.
69integer :: w
70!nb de points récupérés par le calcul de Laskar (e.g. 1e6 années avec un point
71!tous les 100 ans => 10001 points)
72integer,parameter :: nb_points_laskar=5001
73!Les variables TEMPS, INSOL_dec, INSOL_jan, INSOL_fev correspondent au points
74!récupérés directement dans le fichier txt du site de Laskar. INSOL est la
75!moyenne des trois variables précédentes  TEMPS_mod est une échelle de temps
76!recalibrée pour débuter à 0 et être échelonnée en années.
77double precision,dimension(nb_points_laskar) :: TEMPS,TEMPS_mod
78double precision,dimension(nb_points_laskar) :: INSOL,INSOL_dec,INSOL_jan,INSOL_fev
79!yp1_ins et ypn_ins sont les dérivés premières au 1er et dernier point de INSOL.
80double precision :: yp1_ins,ypn_ins
81!var_bidon1 et 2 sont des variables inutiles qui servent à ne pas multiplier le
82!nombre de tableaux contenant les différents points de TEMPS.
83real :: var_bidon1,var_bidon2
84
85real,dimension(nx,ny,mois,ctr,gtr) :: Tm_time_fin_1     ! temperature mensuelle au temps time (non corrige altitude)
86real,dimension(nx,ny,mois,ctr,gtr) :: Pm_time_fin_1     ! precipitation mensuelle au temps time
87
88!interpolation sur surface glace
89!real,allocatable :: nb_points_max(:)
90real,dimension(nx,ny) :: prox_max ! somme des distances pour chaque point
91real,allocatable :: maillage_prox(:,:) ! distance des points autour du point etudie
92real,allocatable :: Hlarge(:,:) ! H avec grille elargie
93integer :: error
94real :: distsq
95real :: real_rayon_ice
96integer :: grid_int ! rayon de recherche en nbr de points de grille
97!integer :: indice
98real,dimension(nx,ny) :: Height_new,Height_old
99real,dimension(nx,ny) :: coeff2_stored
100!real,dimension(nx,ny) :: coeff1_stored,coeff2_stored
101real,dimension(nx,ny,mois,ctr) :: Tm_time_fin_2         ! temperature mensuelle au temps time (non corrige altitude)
102real,dimension(nx,ny,mois,ctr) :: Pm_time_fin_2         ! precipitation mensuelle au temps time
103
104
105!interpolation sur le CO2
106double precision,allocatable :: CO2_time(:), CO2_val(:) ! depend du nbr de points du fichier CO2
107real,dimension(nx,ny,mois) :: Tm_time_fin_3             ! temperature mensuelle au temps time (non corrige altitude)
108real,dimension(nx,ny,mois) :: Pm_time_fin_3             ! precipitation mensuelle au temps time
109
110
111real,dimension(nx,ny,mois,ctr,gtr) :: Tm_surf_mod       ! correction altitude
112real,dimension(nx,ny,mois,ctr,gtr) :: Pm_surf_mod       ! correction altitude
113
114real,dimension(nx,ny,mois) :: Tm_surf                 ! surface temperature (after topo. correction)
115                                                     
116real,dimension(nx,ny,mois) :: Pm_surf                 ! surface precipitation (after topo. correction)
117
118real,dimension(nx,ny,mois) :: lapserate                 ! lapse rate
119real :: psolid=2.                                      ! temp limit between liquid and solid precip
120
121character(len=150)  :: filin                           ! nom temporaire
122character(len=100) :: file_temporel                    ! forcage temporel
123character(len=100),dimension(ctr,gtr,ntr) :: filtr_t   !snapshot temp file name
124character(len=100),dimension(ctr,gtr,ntr) :: filtr_p   ! precip file
125character(len=100),dimension(ctr,gtr) :: surf_ice      ! topo file
126character(len=100),dimension(3) :: orb_file            ! orbite file
127character(len=100) :: co2_file                         ! co2 file
128real,dimension(nx,ny,ctr,gtr) :: Ssnap                 ! topo reference
129
130real,dimension(nx,ny) :: ZS                             !< surface topography above sea level
131
132! specifique greeneem15 : calcul points englace sur Groenland uniquement :
133logical, dimension(nx,ny,2) :: mask_cal                ! masque calotte Groenland (comme dans greeneem/output)
134integer,dimension(2) :: INP                            ! nbr de points englacés
135
136real :: mincoefbmelt                                    ! butoirs pour coefbmshelf
137real :: maxcoefbmelt 
138
139
140
141contains
142
143! 2=lecture des inputs
144!--------------------
145
146subroutine input_clim                                   ! routine qui permet d'initialiser les variables climatiques
147! variables locales
148!-------------------
149
150  implicit none
151 
152  integer :: intr
153  integer :: igtr
154  integer :: ictr
155  integer :: d1,d2
156  integer :: i,j
157!  character(len=100) :: file_ncdf      !< fichier netcdf
158  real*8, dimension(:,:,:), pointer :: data_3D => null() ! donnees lues dans le netcdf
159  real*8, dimension(:,:),pointer :: data_2D => null()    ! donnees lues dans le netcdf
160
161! lecture des fichiers snapshots pour tout geoplace
162! -------------------------------------------------
163write(6,*) 'fichiers snapshots'
164do ictr=1,ctr  ! co2
165   do igtr=1,gtr  ! calotte
166      do intr=1,ntr ! orbitaux
167
168      !temperature
169      filin=trim(dirnameinp)//'forcing/'//trim(filtr_t(ictr,igtr,intr))
170      call Read_ncdf_var('t2m',trim(filin),data_3D)    ! Temperature
171      Tm_fin(:,:,:,ictr,igtr,intr)=data_3D(:,:,:) 
172!      write(6,*) trim(filin)
173!      open(20,file=trim(filin))
174!      do j=1,ny
175!         do i=1,nx
176      !    do J=ny,1,-1
177      !       do I=1,nx
178!            read(20,*) ti, tj, (Tm_fin(i,j,mo,ictr,igtr,intr),mo=1,12)
179!         end do
180!      end do
181!      close(20)
182
183      !precipitation
184      filin=trim(dirnameinp)//'forcing/'//trim(filtr_p(ictr,igtr,intr))
185      call Read_ncdf_var('precip',trim(filin),data_3D)    ! precipitation
186      Pm_fin(:,:,:,ictr,igtr,intr)=data_3D(:,:,:) 
187!      write(6,*) trim(filin)
188!      open(20,file=trim(filin))
189!      do j=1,ny
190!         do i=1,nx
191         !    do J=ny,1,-1
192         !       do I=1,nx
193!            read(20,*) ti, tj, (Pm_fin(i,j,mo,ictr,igtr,intr),mo=1,12)
194!         end do
195!      end do
196!      close(20)
197      end do
198
199       ! topo
200      filin=trim(dirnameinp)//'forcing/'//trim(surf_ice(ictr,igtr))
201      call Read_ncdf_var('TOPO',trim(filin),data_2D)    ! topo
202      Ssnap(:,:,ictr,igtr)=data_2D(:,:) 
203      where(Ssnap(:,:,ictr,igtr).eq.0.0)  ! Pour PLIOMIP niv marin=25m
204         Ssnap(:,:,ictr,igtr)=25.0
205      endwhere
206!      write(6,*) trim(filin)
207!      open(20,file=trim(filin))
208!      do j=1,ny
209!         do i=1,nx
210         !    do J=ny,1,-1
211         !       do I=1,nx
212!            read(20,*) Ssnap(i,j,ictr,igtr)
213!         end do
214!      end do
215!      close(20)
216   end do
217end do
218
219
220!Lecture des parametres orbitaux
221select case(trim(orbital_interpolation))
222case('La2004')
223   
224   print *,'Laskar 2004 orbital parameters'
225   open(70,file=trim(dirnameinp)//'forcing/'//trim(orb_file(1)))
226      do w=1,nb_points_laskar
227         read(70,*) TEMPS(w),INSOL_dec(w)
228      end do
229   close(70)
230   
231   open(71,file=trim(dirnameinp)//'forcing/'//trim(orb_file(2)))
232      do w=1,nb_points_laskar
233         read(71,*) var_bidon1,INSOL_jan(w)
234      end do
235   close(71)
236
237   open(72,file=trim(dirnameinp)//'forcing/'//trim(orb_file(3)))
238      do w=1,nb_points_laskar
239         read(72,*) var_bidon2,INSOL_fev(w)
240      end do
241   close(72)
242
243   do w=1,nb_points_laskar
244      TEMPS_mod(w) = 1000*TEMPS(w)    !+34000000
245      INSOL(w)=(INSOL_dec(w)+INSOL_jan(w)+INSOL_fev(w))/3
246   end do
247
248   !derivée 1er et dernier point
249   yp1_ins=(INSOL(2)-INSOL(1))/(TEMPS_mod(2)-TEMPS_mod(1))
250   ypn_ins=(INSOL(nb_points_laskar)-INSOL(nb_points_laskar-1))/(TEMPS_mod(nb_points_laskar)-TEMPS_mod(nb_points_laskar-1))
251 
252case('constant')
253   print *,'Constant orbital parameters'
254
255case('sinusoid')
256   print *,'Sinusoidal variation of orbital parameters'
257
258end select
259
260!Calcul des poids relatifs lies a la position d'un point par rapport a un autre
261select case(trim(ice_int))
262case('no')
263   print *,"Ice interpolation not prescribed"
264
265case('yes')
266   print *,"Ice interpolation prescribed"
267   
268   select case(trim(interpolation_ice))
269   case('linear')
270      print *,"Linear ice mode"
271
272   case('nonlinear')
273      print *,"Non linear ice mode"
274         
275      ! Calculation of the radius of search in terms of indexes.
276      grid_int=nint(rayon_interpolation_ice/(dx/1000.))
277      real_rayon_ice=real(rayon_interpolation_ice)
278      print *,grid_int,ntr,gtr,ctr,intr,igtr,ictr
279      print*,'dx=',dx
280      print*,'rayon_interpolation_ice=',rayon_interpolation_ice
281      print*,'grid_int=',grid_int
282!      indice=grid_int
283!      do j=1,ny
284!         do i=1,nx
285            !Calculation of the array indice, which stores the number of points
286            !taken into account in the interpolation. Practically overAntarctica,
287            !it is always equal to the number of points within the circle defined
288            !by grid_int, but closer to the edges of the grid, this number is smaller.
289!cdc            indice(i,j)=min(min(i,j),min(nx-(i-1),ny-(j-1)))-1
290!cdc            if (indice(i,j).gt.(grid_int)) then
291!               indice(i,j)=grid_int
292!cdc            endif
293           
294!            if ((j.eq.198).and.(i.eq.25)) then
295!               print *,'time = ',time,'ictr = ',ictr,'np = ',np
296!               print *,'i = ',i,'j = ',j,'indice = ',indice(i,j)
297!            endif
298!            if ((j.eq.197).and.(i.eq.25)) then
299!               print *,'time = ',time,'ictr = ',ictr,'np = ',np
300!               print *,'i = ',i,'j = ',j,'indice = ',indice(i,j)
301!            endif
302!         enddo
303!      enddo
304
305      !Allocation of the array nb_points_max, which stores the maximum number of
306      !points potentially used in the interpolation for every indice value.
307      !allocate(nb_points_max(grid_int),stat=error)
308      !if (error.ne.0) then
309      !   print *,'error: could not allocate memory for array'
310      !   stop
311      !endif
312      !Allocation of the array prox_max. Same as above but for the proximity
313      !parameter.
314!cdc      allocate(prox_max(grid_int),stat=error)
315!      if (error.ne.0) then
316!         print *,'error: could not allocate memory for array'
317!         stop
318!      endif
319      !Allocation of the array maillage_prox, which stores, for every indice
320      !value, the proximity parameter associated with each point within the
321      !radius of search.
322      allocate(maillage_prox(-grid_int:grid_int,-grid_int:grid_int),stat=error)
323      if (error.ne.0) then
324         print *,'error: could not allocate memory for array'
325         stop
326      endif
327!      print*,'Hlarge',grid_int, 1-grid_int, nx+grid_int,ny+grid_int
328
329      allocate(Hlarge(1-grid_int:nx+grid_int,1-grid_int:ny+grid_int),stat=error)
330      if (error.ne.0) then
331         print *,'error: could not allocate memory for array'
332         stop
333      endif
334      Hlarge=0.
335
336      !Initialization at 0 for the three arrays allocated before.
337      !nb_points_max=0
338      prox_max=0
339      maillage_prox=0
340     
341      !Filling of the three arrays. The proximity parameter depends on a ration
342      !of exponential functions.
343!cdc      do k=1,grid_int
344      do d2=-grid_int,grid_int
345        do d1=-grid_int,grid_int
346          distsq=(abs(d1)*(dx/1000.))**2+(abs(d2)*(dx/1000.))**2
347               !print *,'distsq = ',distsq,'real_rayon_ice =',real_rayon_ice
348          if ((sqrt(distsq).le.real_rayon_ice).and.((d1/=0).or.(d2/=0))) then
349                  !nb_points_max(k)=nb_points_max(k)+1
350!cdc                  prox_max(k)=prox_max(k)+100*exp(-((-1*log(0.001))/real_rayon_ice)*(sqrt(distsq)))/exp(-((-1*log(0.001))/real_rayon_ice)*(dx/1000.))
351                  !Indexes displacement to match elements
352              maillage_prox(d1,d2)=100*exp(-((-1*log(0.001))/real_rayon_ice)*(sqrt(distsq)))/exp(-((-1*log(0.001))/real_rayon_ice)*(dx/1000.))
353                  !print *,'nb_points_max = ',nb_points_max(k),'prox_max = ',prox_max(k),'maillage_prox = ',maillage_prox(k,d1+(k+1),d2+(k+1))
354          endif
355        enddo
356      enddo
357!cdc      enddo
358
359     
360      prox_max(:,:)=SUM(maillage_prox(:,:))
361!      PRINT *,'maillage_prox =',maillage_prox(:,:)
362!      PRINT*,'sum(maillage_prox)',SUM(maillage_prox(:,:))
363!      PRINT*,'fin calcul maillage prox'
364   end select
365end select
366
367Height_new=0.
368Height_old=0.
369!coeff1_stored=0.
370coeff2_stored=0.
371
372select case(trim(CO2_function))
373case('constant')
374   print *,'CO2 kept constant at ',CO2_value
375
376case('variable')
377   print *,'The CO2 will vary with time'
378   open(75,file=trim(dirnameinp)//'forcing/'//trim(co2_file))
379   read(75,*) nb_slope_dis_CO2
380   allocate(CO2_time(nb_slope_dis_CO2),stat=error)
381   if (error.ne.0) then
382      print *,'error: could not allocate memory for array CO2_time'
383      stop
384   endif
385   allocate(CO2_val(nb_slope_dis_CO2),stat=error)
386   if (error.ne.0) then
387      print *,'error: could not allocate memory for array CO2_val'
388      stop
389   endif
390   do w=1,nb_slope_dis_CO2
391      read(75,*) CO2_time(w),CO2_val(w)
392   end do
393   close(75)
394   !print *,'CO2_time = ',CO2_time(:),'CO2_val = ',CO2_val(:)
395end select
396
397
398! initilisation region Groenland pour calcul nbr points englaces
399mask_cal(:,:,:)=.false.
400mask_cal(:,:,1)=.true.
401do j=1,ny
402   do i=1,nx
403! -- Groenland
404      IF ( (((xlong(i,j).ge.290).AND.(xlong(i,j).le.350)).AND. &
405           ((ylat(i,j).ge.79.5).AND.(ylat(i,j).le.85))).OR.      &
406           (((xlong(i,j).ge.286).AND.(xlong(i,j).lt.345)).AND. &
407           ((ylat(i,j).ge.75).AND.(ylat(i,j).le.79.5))).OR.      &
408           (((xlong(i,j).ge.300).AND.(xlong(i,j).le.345)).AND. &
409           ((ylat(i,j).ge.67).AND.(ylat(i,j).le.75))).OR.    &
410           (((xlong(i,j).ge.305).AND.(xlong(i,j).le.330)).AND. &
411           ((ylat(i,j).ge.55).AND.(ylat(i,j).le.67))) ) mask_cal(i,j,2)=.true.
412   enddo
413enddo
414
415
416end subroutine input_clim
417!--------------------------------------------------------------------------------
418!subroutine input_climat_ref
419! quand on traite en absolu, pas besoin du climat de reference
420
421!end subroutine input_climat_ref
422
423
424
425SUBROUTINE init_forclim
426
427! fichiers snapshots
428  namelist/snap_forcage_mois_JB/filtr_t,filtr_p,Seuil_haut,Seuil_bas,summorb,palier_ice,surf_ice,palier_CO2,orb_file,co2_file  ! ce bloc est a dupliquer pour chaque snapshot en changeant                                                     ! la numerotation. ntr snapshots
429
430! forcage temporel
431!------------------
432  namelist/forc_temporel/file_temporel,mincoefbmelt,maxcoefbmelt
433
434! lecture par namelist
435!---------------------
436! formats pour les ecritures dans 42
437428 format(A)
438  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
439  read(num_param,snap_forcage_mois_JB)
440  write(num_rep_42,428)'!___________________________________________________________' 
441  write(num_rep_42,428) '&snap_forcage_mois_JB                            ! module climat-forcage-insolation_mod'
442  write(num_rep_42,'(A,A)')   'filtr_t = ', filtr_t
443  write(num_rep_42,'(A,A)')   'filtr_p = ', filtr_p
444  write(num_rep_42,'(A,2(f7.1,","))')   'Seuil_haut = ', Seuil_haut(:)
445  write(num_rep_42,'(A,2(f7.1,","))')   'Seuil_bas = ', Seuil_bas(:)
446  write(num_rep_42,'(A,2(f5.1,","))')   'summorb = ', summorb(:)
447  write(num_rep_42,'(A,2(f7.1,","))')   'palier_ice = ', palier_ice(:,:)
448  write(num_rep_42,'(A,A)')   'surf_ice = ', surf_ice
449  write(num_rep_42,'(A,2(f6.1,","))')   'palier_CO2 = ', palier_CO2(:)
450  write(num_rep_42,'(A,A)')   'orb_file = ', orb_file
451  write(num_rep_42,'(A,A)')   'co2_file = ', co2_file
452  write(num_rep_42,*)'/'                     
453  write(num_rep_42,428) '! fichiers temperature et precip : 12 mois ordre : '
454  write(num_rep_42,428) '! 1/ CO2, 2/ regrowth, 3/ insol. insolmax_regrowthmaw_co2max --> insolmin_noice_co2min '
455  write(num_rep_42,*)
456
457!  do i=1,ntr
458! glaciaire
459!     write(filtr_t(i),'(A,i3,A)') filtr_t1(1:32),int(ttr(i)),filtr_t1(36:50)
460!     write(filtr_p(i),'(A,i3,A)') filtr_p1(1:34),int(ttr(i)),filtr_p1(38:52)
461!     write(filtr_t(i),'(A)') filtr_t
462!     write(filtr_p(i),'(A)') filtr_p
463!  enddo
464
465
466call lect_lapserate_months                               ! lit les lasperate mensuels
467                                                         ! pour une version spatialisee ecrire une autre routine
468! fichiers donnant l'evolution temporelle
469! ---------------------------- ------------
470
471rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
472read(num_param,forc_temporel)
473
474write(num_rep_42,428)'!___________________________________________________________' 
475write(num_rep_42,428) '&forc_temporel                            ! module climat_forcage_mois_mod'
476write(num_rep_42,'(A,A)')   'file_temporel =', file_temporel
477write(num_rep_42,*)     'mincoefbmelt      =', mincoefbmelt 
478write(num_rep_42,*)     'maxcoefbmelt      =', maxcoefbmelt 
479write(num_rep_42,*)'/'                     
480write(num_rep_42,428) '!fichier forcage temporel pour snapshot'
481write(num_rep_42,*)
482
483
484end subroutine init_forclim
485!---------------------------------------------------------------------
486
487!forcage climatique au cours du temps
488
489subroutine forclim
490
491implicit none
492
493!integer l              ! dumm index for loops on snapshots files  l=ITR,NTR-1
494!cdc integer itr            ! index of the current snapshot file (change with time)
495
496integer mo
497integer :: ictr
498integer :: igtr
499
500
501!Définition des variables propres au calcul de la fonction d'Insolation
502!y2_ins est le tableau contenant l'ensemble des dérivées secondes au points considérés (cf calcul de
503!l'interpolation par spline cubique).
504double precision,dimension(nb_points_laskar) :: y2_ins
505!insolval est la valeur de l'insolation d'été au pas de temps considere
506double precision :: insolval
507
508!initialisation des paramètres pour l'interpolation entre les fichiers de glace
509!real,allocatable :: sub_H(:,:)
510!real,allocatable :: sub_maillage_prox(:,:)
511!real :: nb_points,prox
512real :: prox
513!real,dimension(nx,ny) :: coeff1,coeff2
514real,dimension(nx,ny) :: coeff2
515integer :: m1,m2
516integer,dimension(nx,ny) :: Mask,Mask2
517integer :: update_H
518integer :: i,j,l
519
520!Definition des variables pour le calcul de la fonction de CO2
521!double precision :: coefCO2
522
523!*****************************************************************************
524! triple interpolation orbite, glace, CO2 + correction altitude
525!*****************************************************************************
526
527
528select case(trim(CO2_function))
529case('constant')
530   CO2SnapMin=1
531   CO2SnapMax=1
532
533case('variable')
534   if (time.le.CO2_time(1)) then
535      CO2_value=CO2_val(1)
536   else if (time.ge.CO2_time(nb_slope_dis_CO2)) then
537      CO2_value=CO2_val(nb_slope_dis_CO2)
538   else
539      l=1
540      do while (time.gt.CO2_time(l))
541         l = l+1
542      end do
543     
544      CO2_value=(((CO2_val(l)-CO2_val(l-1))/(CO2_time(l)-CO2_time(l-1)))*(time-CO2_time(l-1))+CO2_val(l-1))
545 
546! pour ne pas depasser les valeurs CO2 des 2 snapshots min et max CO2
547      if (CO2_value.lt.palier_CO2(ctr)) then
548         CO2_value=palier_CO2(ctr)
549      else if (CO2_value.gt.palier_CO2(1)) then
550         CO2_value=palier_CO2(1)
551      endif
552
553   endif
554
555   l=1
556   do while (CO2_value.lt.palier_CO2(l+1))
557      l=l+1
558   enddo
559   CO2SnapMin=l
560   CO2SnapMax=l+1
561   print *,'time = ',time,'CO2_value = ',CO2_value, 'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax
562
563!   if (CO2_value.ge.360) then
564!      CO2SnapMin=1
565!      CO2SnapMax=2
566!   else if ((coefCO2.lt.360).and.(coefCO2.ge.280)) then
567!      CO2SnapMin=2
568!      CO2SnapMax=3
569!   else if (coefCO2.lt.280) then
570!      CO2SnapMin=3
571!      CO2SnapMax=4
572!   endif
573
574   !print *,'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax
575end select
576
577!*************************************************************************
578!*     triple interpolation orbite, glace, CO2 + correction altitude     *
579!*************************************************************************
580
581
582!Nouvelle fonction d'interpolation des param orbitaux qui permet de récupérer
583!l'insolation réelle de Laskar 2004
584
585
586select case(trim(orbital_interpolation))
587case('La2004')
588   
589   call spline(TEMPS_mod,INSOL,nb_points_laskar,yp1_ins,ypn_ins,y2_ins)
590
591   call splint(TEMPS_mod,INSOL,y2_ins,nb_points_laskar,insolval)
592 
593! Dans massud_param_list.dat, les snapshots doivent etre dans l'ordre
594! decroissant partant du plus chaud
595   if (insolval.gt.summorb(1)) then
596      Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,1)
597      Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,1)
598   else if (insolval.lt.summorb(ntr)) then
599      Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,ntr)
600      Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,ntr)
601   else
602      l=1
603      do while (insolval.lt.summorb(l))
604         l = l+1
605      enddo
606      factsin=(insolval-summorb(l))/(summorb(l-1)-summorb(l))
607      Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,l-1)*factsin+Tm_fin(:,:,:,:,:,l)*(1-factsin)
608      Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,l-1)*factsin+Pm_fin(:,:,:,:,:,l)*(1-factsin)
609   endif
610
611
612case('constant')
613   !Tm_time_fin_1(:,:,:,:,:)=(Tm_fin(:,:,:,:,:,1)+Tm_fin(:,:,:,:,:,2))/2
614   !Pm_time_fin_1(:,:,:,:,:)=(Pm_fin(:,:,:,:,:,1)+Pm_fin(:,:,:,:,:,2))/2
615   Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,1)
616   Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,1)
617
618
619case('sinusoid')
620   orbite=(1+sin(2*pi*time/40000.))/2.
621   Tm_time_fin_1(:,:,:,:,:)=orbite*Tm_fin(:,:,:,:,:,1)+(1-orbite)*Tm_fin(:,:,:,:,:,2)
622   Pm_time_fin_1(:,:,:,:,:)=orbite*Pm_fin(:,:,:,:,:,1)+(1-orbite)*Pm_fin(:,:,:,:,:,2)
623
624end select
625
626
627!do ictr=CO2SnapMin,CO2SnapMax,1
628!   do igtr=1,gtr
629!      print *,surf_ice(ictr,igtr)
630!   enddo
631!enddo
632
633
634!Correction d'altitude pour les états
635!il faut que Ssnap soit un tableau avec les différentes hauteurs de glace selon la
636!simulation (no ice, med ice, full ice et 2x, 2.5x, 3x)
637do j=1,ny
638   do i=1,nx
639      Zs(i,j)=max(sealevel_2d(i,j),S(i,j))
640      do ictr=CO2SnapMin,CO2SnapMax,1
641         do igtr=1,gtr
642            if (Zs(i,j).ge.Ssnap(i,j,ictr,igtr)) then
643               do mo=1,mois
644                  ! correction d'altitude
645                  Tm_surf_mod(i,j,mo,ictr,igtr)=-lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,ictr,igtr)) &
646                                            +Tm_time_fin_1(i,j,mo,ictr,igtr)
647
648                  Pm_surf_mod(i,j,mo,ictr,igtr)=Pm_time_fin_1(i,j,mo,ictr,igtr)*exp(0.05*(Tm_surf_mod(i,j,mo,ictr,igtr) &
649                                            -Tm_time_fin_1(i,j,mo,ictr,igtr)))
650
651               end do
652            else if (Zs(i,j).lt.Ssnap(i,j,ictr,igtr)) then 
653               do mo=1,mois
654                  ! correction d'altitude avec condition pour que T ne soit
655                  ! pas supérieure à la T de l'état sans glace corrigé
656                  Tm_surf_mod(i,j,mo,ictr,igtr)=min(-lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,ictr,igtr))+Tm_time_fin_1(i,j,mo,ictr,igtr), &
657                                            -lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,ictr,gtr))+Tm_time_fin_1(i,j,mo,ictr,gtr))
658
659                  Pm_surf_mod(i,j,mo,ictr,igtr)=Pm_time_fin_1(i,j,mo,ictr,igtr)*exp(0.05*(Tm_surf_mod(i,j,mo,ictr,igtr) &
660                                            -Tm_time_fin_1(i,j,mo,ictr,igtr)))
661               end do
662            endif
663         end do
664      end do
665   end do
666end do
667!print*,'debug Tm_surf_mod corr vert ', Tm_surf_mod(25,198,1,:,:)
668
669! calcul du nbr de points glace pose :
670do l=1,2
671   INP(l) = count(mask_cal(:,:,l).and.(H(:,:).gt.2.).and..not.flot(:,:))
672enddo
673np = INP(2)
674! np = count((H(:,:).gt.2.).and..not.flot(:,:)) ! version std
675!Interpolation glace
676surf_test=max(0,np)
677Height_new=H
678print *,'time = ',time
679!print *,'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax
680!print *,'Seuil_haut(CO2SnapMin) = ',Seuil_haut(CO2SnapMin),'Seuil_haut(CO2SnapMax) = ',Seuil_haut(CO2SnapMax)
681!print *,'Seuil_bas(CO2SnapMin) = ',Seuil_bas(CO2SnapMin),'Seuil_bas(CO2SnapMax) = ',Seuil_bas(CO2SnapMax)
682print *,'surf_test = ',surf_test
683!print *,'coeff1_stored(154,154) = ',coeff1_stored(154,154),'coeff2_stored(154,154) = ',coeff2_stored(154,154)
684!print *,'coeff1_stored(184,184) = ',coeff1_stored(184,184),'coeff2_stored(184,184) = ',coeff2_stored(184,184)
685!print *,'coeff1_stored(200,117) = ',coeff1_stored(200,117),'coeff2_stored(200,117) = ',coeff2_stored(200,117)
686
687select case(trim(ice_int))
688case('no')
689   Tm_time_fin_2(:,:,:,:)=Tm_surf_mod(:,:,:,:,1)
690   Pm_time_fin_2(:,:,:,:)=Pm_surf_mod(:,:,:,:,1)
691
692case('yes')
693   do ictr=CO2SnapMin,CO2SnapMax,1
694      if (surf_test.ge.Seuil_haut(ictr)) then
695         Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,1)
696         Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,1)
697
698         update_H=1
699      else if (surf_test.le.Seuil_bas(ictr)) then
700         Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,gtr)
701         Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,gtr)
702         
703         update_H=0
704      else
705         l=1
706         do while (surf_test.lt.palier_ice(ictr,l))
707            l = l+1
708         end do
709     
710         select case(trim(interpolation_ice))
711         case('linear')
712            retroicecoeff=(surf_test-palier_ice(ictr,l))/((palier_ice(ictr,l-1)-palier_ice(ictr,l)))
713            Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,l-1)*retroicecoeff+Tm_surf_mod(:,:,:,ictr,l)*(1-retroicecoeff)
714            Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,l-1)*retroicecoeff+Pm_surf_mod(:,:,:,ictr,l)*(1-retroicecoeff)
715
716         case('nonlinear')
717            !Important: the distinction between ice and vegetation points used
718            !before has been removed because it would end up with giving the
719            !same interpolation. In this code, the calculation is done only over
720            !surrounding ice points.
721            where (((Height_new > 5.).and.(Height_old < 5.)).or.((Height_new < 5.).and.(Height_old > 5.)))
722               Mask = 1.
723            elsewhere
724               Mask = 0.
725            end where
726           
727            !print *,'Mask = ',Mask
728
729!            allocate(sub_maillage_prox(2*indice+1,2*indice+1),stat=error)
730!            if (error.ne.0) then
731!               print *,'error: could not allocate memory for array'
732!            endif
733
734            Mask2 = 0.
735            do j=1,ny
736               do i=1,nx
737                  if (Mask(i,j).eq.1) then
738                     !print *,'i = ',i,'j = ',j,'Mask = ',Mask(i,j)
739                     do m2=max(1,j-grid_int),min(j+grid_int,ny)  !cdc
740                        do m1=max(i-grid_int,1),min(i+grid_int,nx)  !cdc
741                           Mask2(m1,m2) = 1.
742                        end do
743                     end do
744                  endif
745               end do
746            end do
747           
748            !print *,'Mask2 = ',Mask2
749            Hlarge(1:nx,1:ny)=H(:,:)
750            do j=1,ny
751               do i=1,nx
752                  if (Mask2(i,j).eq.1) then
753                     !The number of ice points nb_points is initialized at 0
754                     !Same for the proximity parameter, prox.
755                     !nb_points=0
756!                     prox=0
757                     !Case where indice = 0, on the edges of the grid.
758                     !We actually dont care here because Antarctica is far from the
759                     !edges of the grid.
760                 
761                     !if ((j.eq.154).and.(i.eq.154)) then
762                     !   print *,'indice =',indice(i,j)
763                     !endif
764
765                     !If indice = 0, no interpolation, the climate used comes from
766                     !the deglaciated simulation (practically, we dont care about
767                     !this because it concerns only the 1st and last rows and
768                     !columns of the grid, so not Antarctica.... Could be
769                     !problematic in other contexts)
770                     !if (indice(i,j).eq.0) then
771                     !   Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)
772                     !   Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)
773                       
774                     !General case
775                     !else   
776                        !Allocation of the array sub_H, which will be used to count
777                        !the number of ice points around the focused point.
778                        !allocate(sub_H(2*indice(i,j)+1,2*indice(i,j)+1),stat=error)
779                        !if (error.ne.0) then
780                        !   print *,'error: could not allocate memory for array'
781                        !endif
782                        !sub_H is initialized at 1 everywhere except on the
783                        !focused point which is not taken into account.
784                        !sub_H=1
785                        !sub_H(indice(i,j)+1,indice(i,j)+1)=0
786                     
787                        !where (maillage_prox(indice(i,j),:,:) == 0.)
788                        !   sub_H(:,:)=0.
789                        !end where
790
791                        !Calculation of the number of ice points present around the
792                        !focused one.
793                        !nb_points=sum(sub_H,mask=H(i-indice(i,j):i+indice(i,j):1,j-indice(i,j):j+indice(i,j):1).ge.5)
794                     
795                        !if ((j.eq.117).and.(i.eq.200)) then
796                        !   print *,'time =',time,'sub_H =',sub_H(:,:),'nb_points =',nb_points,'H =',H(i-indice(i,j):i+indice(i,j):1,j-indice(i,j):j+indice(i,j):1)
797                        !   print *,'time =',time,'nb_points =',nb_points
798                        !endif
799                     
800                        !deallocate(sub_H)
801
802                        !Allocation of the array sub_maillage_prox, which is
803                        !extracted from the array maillage_prox (calculated above
804                        !in input_clim) accordingly to the position of the focused
805                        !point (ie the value of indice).
806!                        allocate(sub_maillage_prox(2*indice(i,j)+1,2*indice(i,j)+1),stat=error)
807!                        if (error.ne.0) then
808!                           print *,'error: could not allocate memory for array'
809!                        endif
810                        !sub_maillage_prox is initialized at 0 everywhere.
811!                        sub_maillage_prox=0
812                     
813                        !sub_maillage_prox is then filled with the corresponding
814                        !values from maillage_prox.
815!                        do k2=1,2*indice+1
816!                           do k1=1,2*indice+1
817!                              sub_maillage_prox(k1,k2)=maillage_prox(k1,k2)
818!                           end do
819!                        end do
820                        prox=0.
821                        !Calculation of the proximity parameter.
822!                        do k2=max(j-indice,1),min(j+indice,ny)
823!                           do k1=max(i-indice,1),min(i+indice,nx)
824!                              if (H(k1,k2).ge.5) prox=prox+sub_maillage_prox(k1,k2)
825!                           enddo
826!                        enddo
827
828!                        maillage_prox_ij(i-grid_int:i+grid_int:1,j-grid_int:j+grid_int:1)=maillage_prox
829
830
831
832                        prox=sum(maillage_prox,mask=Hlarge(i-grid_int:i+grid_int:1,j-grid_int:j+grid_int:1).ge.5)
833                       
834!                        if ((j.eq.198).and.(i.eq.25)) then
835!                           print *,'maillage_prox =',maillage_prox(:,:),'prox =',prox
836!                           print *,'prox_max =',prox_max(i,j)
837!                        endif
838   
839!                        deallocate(sub_maillage_prox)
840                     
841                        !Calculation of the relative coefficient (<1) for the
842                        !number of points and the proximity.
843                        !coeff1(i,j)=nb_points/nb_points_max(indice(i,j))
844                        coeff2(i,j)=prox/prox_max(i,j)
845!                        PRINT*,'i=',i,'j=',j
846!                        PRINT*,'prox=',prox
847!                        PRINT*,'prox_max=',prox_max(i,j)
848!                        PRINT*,'coeff2=',coeff2(i,j)
849!                        DO m1=j+grid_int,j-grid_int,-1
850!                          PRINT*,'Hlarge',Hlarge(i-grid_int:i+grid_int:1,m1)
851!                        ENDDO
852!                        PRINT*
853
854                        !if (((j.eq.154).and.(i.eq.154)).or.((j.eq.184).and.(i.eq.184)).or.((j.eq.117).and.(i.eq.200))) then
855                        !       print *,'time = ',time,'i = ',i,'j = ',j,'coeff1 = ',coeff1(i,j),'coeff2 = ',coeff2(i,j)
856                        !endif
857
858                        !Interpolation:
859                        !Tm_surf_mod_noice*((1-a*coeff1)+(1-b*coeff2))/2+Tm_surf_mod_ice*(a*coeff1+b*coeff2)/2
860                        !Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)*((1-(2*poids_nb_points)*coeff1(i,j))+(1-(2*poids_proximity)*coeff2(i,j)))/2+Tm_surf_mod(i,j,:,ictr,l-1)*((2*poids_nb_points)*coeff1(i,j)+(2*poids_proximity)*coeff2(i,j))/2
861                        !Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)*((1-(2*poids_nb_points)*coeff1(i,j))+(1-(2*poids_proximity)*coeff2(i,j)))/2+Pm_surf_mod(i,j,:,ictr,l-1)*((2*poids_nb_points)*coeff1(i,j)+(2*poids_proximity)*coeff2(i,j))/2
862                       
863                        !Pour après
864                        Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)*(1-coeff2(i,j))+Tm_surf_mod(i,j,:,ictr,l-1)*coeff2(i,j)
865                        Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)*(1-coeff2(i,j))+Pm_surf_mod(i,j,:,ictr,l-1)*coeff2(i,j)
866                       
867                        !coeff1_stored(i,j)=coeff1(i,j)
868                        coeff2_stored(i,j)=coeff2(i,j)
869                        !print *,'coeff1_stored(154,154) = ',coeff1_stored(154,154),'coeff2_stored(154,154) = ',coeff2_stored(154,154)
870                        !print *,'coeff1_stored(184,184) = ',coeff1_stored(184,184),'coeff2_stored(184,184) = ',coeff2_stored(184,184)
871                        !print *,'coeff1_stored(200,117) = ',coeff1_stored(200,117),'coeff2_stored(200,117) = ',coeff2_stored(200,117)
872                     !endif
873                  else
874                     !Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)*((1-(2*poids_nb_points)*coeff1_stored(i,j))+(1-(2*poids_proximity)*coeff2_stored(i,j)))/2+Tm_surf_mod(i,j,:,ictr,l-1)*((2*poids_nb_points)*coeff1_stored(i,j)+(2*poids_proximity)*coeff2_stored(i,j))/2
875                     !Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)*((1-(2*poids_nb_points)*coeff1_stored(i,j))+(1-(2*poids_proximity)*coeff2_stored(i,j)))/2+Pm_surf_mod(i,j,:,ictr,l-1)*((2*poids_nb_points)*coeff1_stored(i,j)+(2*poids_proximity)*coeff2_stored(i,j))/2
876                     !Pour après
877                     Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)*(1-coeff2_stored(i,j))+Tm_surf_mod(i,j,:,ictr,l-1)*coeff2_stored(i,j)
878                     Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)*(1-coeff2_stored(i,j))+Pm_surf_mod(i,j,:,ictr,l-1)*coeff2_stored(i,j)
879                  endif
880
881               end do
882            end do
883!            deallocate(sub_maillage_prox)
884         end select
885
886         update_H=1
887      endif
888   end do
889end select
890!print*,'debug Tm_time_fin_2 ', Tm_time_fin_2(25,198,1,:)
891!print*,'debug coeff2_stored ', coeff2_stored(25,198)
892!print*,'debug grid_int ', grid_int
893
894
895!do j=1,ny
896!   do i=1,nx
897!      if (coeff1_stored(i,j).gt.0) then
898!         print *,'i = ',i,'j = ',j,'coeff1_stored = ',coeff1_stored(i,j)
899!      endif
900!      if (coeff2_stored(i,j).gt.0) then
901!         print *,'i = ',i,'j = ',j,'coeff2_stored = ',coeff2_stored(i,j)
902!      endif
903!   end do
904!end do
905
906if (update_H.eq.0) then
907   Height_old=0.
908   print *,'No update of H'
909else if (update_H.eq.1) then
910   Height_old=Height_new
911   print *,'H updated'
912else
913   print *,'Ice sheet topography error. Exit'
914   stop
915endif
916!print *,'time = ',time,'TEMP(200,117) = ',Tm_time_fin_2(200,117,:,ictr),'PRECIP(200,117) = ',Pm_time_fin_2(200,117,:,ictr)
917
918!fonction de décroissance du CO2
919!CO2 max en 1er puis décroissance dans l'ordre
920
921select case(trim(CO2_function))
922case('constant')
923   Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,1)
924   Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,1)
925
926
927case DEFAULT
928   print *,CO2_value
929   if (CO2_value.ge.palier_CO2(1)) then
930      Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,1)
931      Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,1)
932   else if (CO2_value.le.palier_CO2(ctr)) then
933      Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,ctr)
934      Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,ctr)
935   else
936      l=1
937      do while (CO2_value.lt.palier_CO2(l))
938         l = l+1
939      end do
940      select case(trim(interpolation_CO2))
941      case('linear')
942         factCO2=(CO2_value-palier_CO2(l))/((palier_CO2(l-1)-palier_CO2(l)))
943      case('logarithmic')
944         factCO2=log(CO2_value/palier_CO2(l))/log(palier_CO2(l-1)/palier_CO2(l))
945      end select
946      Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,l-1)*factCO2+Tm_time_fin_2(:,:,:,l)*(1-factCO2)
947      Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,l-1)*factCO2+Pm_time_fin_2(:,:,:,l)*(1-factCO2)
948   endif
949end select
950!print*,'debug Tm_time_fin_3 ', Tm_time_fin_3(25,198,1)
951!print *,time,h(196,145),uxbar(196,145),uybar(196,145),bm(196,145)
952
953do i=1,nx
954   do j=1,ny
955      do mo=1,mois
956         Tm_surf(i,j,mo)=Tm_time_fin_3(i,j,mo)
957         Pm_surf(i,j,mo)=Pm_time_fin_3(i,j,mo)
958      end do
959   end do
960end do
961!print*,'debug Tm_surf = Tm_3 ', Tm_surf(25,198,1)
962
963!coefbmshelf coefficient pour la fusion basale sous les ice shelves
964coefbmshelf=1.0
965coefbmshelf=max(coefbmshelf,mincoefbmelt)
966coefbmshelf=min(coefbmshelf,maxcoefbmelt)
967
968
969
970
971!*********************************
972!Correction d'altitude à vérifier
973!*********************************
974
975
976!  Correction d'altitude pour la temperature y compris sur les lacs
977! Zs est l'altitude de la surface qu'elle soit mer, glace ou lac
978
979!do j=1,ny
980!   do i=1,nx
981!      ZS(I,J)=max(slv(i,j),S(I,J))
982!      do mo=1,mois
983!         Tm_surf(i,j,mo)= - lapserate(i,j,mo) * (Zs(i,j)-Ssnap(i,j)) & ! correction d'altitude T
984!              + Tm_time(i,j,mo)
985!
986!         Pm_surf(i,j,mo)= Pm_time(i,j,mo)*exp(0.05*(Tm_surf(i,j,mo)-Tm_time(i,j,mo)))
987!
988!      end do
989!   end do
990!end do
991
992
993do mo=1,mois
994   Tmois(:,:,mo)=Tm_surf(:,:,mo)
995enddo
996!print*,'debug Tmois ', Tmois(25,198,1)
997
998! calcul de Tann et Tjuly pour les sorties :
999Tann(:,:)=sum(Tmois,dim=3)/12.     ! moy annuelle
1000Tjuly(:,:)=Tmois(:,:,7)            ! temp juillet
1001!print*,'debug Tann ', Tann(25,198)
1002
1003acc(:,:)=sum(Pm_surf,dim=3,mask=Tmois < psolid) ! /12.
1004acc(:,:)=acc(:,:)*1000./ro
1005
1006! pas beau ATTENTION
1007! pb plantage zone spitzberg bord de grille !!!!!
1008!do j=180,190
1009!   do i=123,126
1010!      acc(i,j)=0.
1011!      print*,'bsoc(i,j)',i,j,bsoc(i,j)
1012!   enddo
1013!enddo
1014
1015END subroutine forclim
1016
1017
1018
1019!************************************************************************
1020! Numerical Recipes
1021! interpolation spline cubique
1022
1023!Récupéré grace à Christophe. Modifié 17.04.13 par JB
1024SUBROUTINE splint(xa,ya,y2a,n,y)
1025  INTEGER,intent(in) :: n
1026  double precision,dimension(n),intent(in) :: xa,y2a,ya
1027  double precision,intent(out) :: y
1028! Calculates the cubic spline interpolation.
1029! Given 2 arrays of dimension n, xa and ya, and y2a, the second derivative
1030! of the function ya at any of the n points, it computes the interpolation for
1031! the array y of dimension nmax.
1032! for example, if n is 10001 (e.g. 1e6 years computed with Laskar algorithm and
1033! a sampling step of 100 years), then nmax would be 1e6 because GRISLI has a
1034! sampling step of 1 year.
1035  INTEGER l
1036  REAL a,b,h
1037! We will find the right place in the table by means of bisection.
1038
1039! do i=1,nmax
1040! print *,x(i)
1041! write (*,*) 'Press Enter to Continue'
1042! read (*,*)
1043! enddo
1044
1045!l=1
1046!do j=1,nmax
1047!if (j==1) then
1048!   l=2
1049!else if (j==nmax) then
1050!   l=n
1051!else if (x(j).gt.xa(l)) then
1052!   l=l+1
1053!endif
1054
1055
1056if (time==0) then
1057   l=2
1058else if (time==xa(n)) then
1059   l=n
1060else
1061   l=1
1062   do while (time.gt.xa(l))
1063   l=l+1
1064   enddo
1065endif
1066
1067
1068h=xa(l)-xa(l-1)
1069a=(xa(l)-time)/h
1070b=(time-xa(l-1))/h
1071y=a*ya(l-1)+b*ya(l)+((a**3-a)*y2a(l-1)+(b**3-b)*y2a(l))*(h**2)/6.
1072
1073! print *,time,l,xa(l-1),xa(l),ya(l-1),ya(l),h,a,b,y
1074! write (*,*) 'Press Enter to Continue'
1075! read (*,*)
1076
1077return
1078
1079END SUBROUTINE splint
1080
1081
1082
1083! Calculates the 2nd derivative of the y function for any points
1084! Recupere grace à Christophe. Modifié par JB 18.04.13
1085SUBROUTINE spline(x,y,n,yp1,ypn,y2)
1086  INTEGER,intent(in) :: n
1087  double precision,dimension(n), intent(in) :: x,y
1088  double precision,dimension(n), intent(out) :: y2
1089  double precision,intent(in) :: yp1,ypn
1090
1091  INTEGER i,k
1092  double precision :: p,qn,sig,un
1093  double precision,dimension(n) :: u
1094
1095if (yp1.gt..99e30) then ! The lower boundary condition is set either to 0
1096   y2(1)=0.
1097   u(1)=0.
1098else ! or else to have a specified first derivative.
1099   y2(1)=-0.5
1100   u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
1101endif
1102do i=2,n-1
1103! This is the decomposition loop of the tridiagonal algorithm. y2 and u are used
1104! for temporary storage of the decomposed factors.
1105
1106   sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
1107   p=sig*y2(i-1)+2.
1108   y2(i)=(sig-1.)/p
1109   u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
1110        /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
1111enddo
1112
1113if (ypn.gt..99e30) then ! The upper boundary condition is set eith
1114   qn=0.
1115   un=0.
1116else ! or else to have a specified first derivative.
1117   qn=0.5
1118   un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
1119endif
1120y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
1121do k=n-1,1,-1 ! This is the backsubstitution loop of the tridiagonal algorithm
1122   y2(k)=y2(k)*y2(k+1)+u(k)                             
1123enddo
1124return
1125
1126END SUBROUTINE spline
1127
1128
1129
1130
1131
1132
1133
1134!------------------------------------------------------------------------------------------------------------
1135subroutine lect_lapserate_months  ! lapserates mensuels mais uniformes spatialement
1136
1137implicit none
1138real,dimension(12) :: lect_lapse       ! pour la lecture
1139integer :: i,j
1140
1141namelist/lapse_month/lect_lapse
1142
1143! lecture de la namelist
1144
1145! formats pour les ecritures dans 42
1146428 format(A)
1147
1148rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
1149read(num_param,lapse_month)
1150
1151write(num_rep_42,428)'!___________________________________________________________' 
1152write(num_rep_42,428) '&lapse_month                                          !  module climat_forcage_mois_mod'
1153write(num_rep_42,'(A,12(f0.2,","))')  'lapse_month  = ', lect_lapse
1154write(num_rep_42,*)'/'                     
1155write(num_rep_42,428) '! laspe rates janvier -> decembre en deg/km'
1156
1157
1158! pour repasser en deg/m et copier dans lapserate
1159
1160do j=1,ny
1161   do i=1,nx
1162      lapserate(i,j,:)=lect_lapse(:)/1000.
1163   end do
1164end do
1165end subroutine lect_lapserate_months
1166!--------------------------------------------------------------------------------------------------------
1167
1168end module climat_forcage_insolation_mod
Note: See TracBrowser for help on using the repository browser.