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

Last change on this file since 19 was 19, checked in by dumas, 9 years ago

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