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

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

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