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

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

climat-forcage-insolation_mod.f90 : Lecture du nombre de ligne du fichier CO2 a partir de la premiere ligne du fichier de forcage (plus de modification a faire en dur dans le code)
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!real :: poids_nb_points=0.5
38!real :: poids_proximity=0.5
39character(len=15) :: CO2_function='variable'   !Choose between variable and constant
40integer :: nb_slope_dis_CO2  ! 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,allocatable :: CO2_time(:), CO2_val(:) ! depend du nbr de points du fichier CO2
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   open(75,file=trim(dirnameinp)//'forcing/'//trim(co2_file))
392   read(75,*) nb_slope_dis_CO2
393   allocate(CO2_time(nb_slope_dis_CO2),stat=error)
394   if (error.ne.0) then
395      print *,'error: could not allocate memory for array CO2_time'
396      stop
397   endif
398   allocate(CO2_val(nb_slope_dis_CO2),stat=error)
399   if (error.ne.0) then
400      print *,'error: could not allocate memory for array CO2_val'
401      stop
402   endif
403   do w=1,nb_slope_dis_CO2
404      read(75,*) CO2_time(w),CO2_val(w)
405   end do
406   close(75)
407   !print *,'CO2_time = ',CO2_time(:),'CO2_val = ',CO2_val(:)
408end select
409
410
411! initilisation region Groenland pour calcul nbr points englaces
412mask_cal(:,:,:)=.false.
413mask_cal(:,:,1)=.true.
414do j=1,ny
415   do i=1,nx
416! -- Groenland
417      IF ( (((xlong(i,j).ge.290).AND.(xlong(i,j).le.350)).AND. &
418           ((ylat(i,j).ge.79.5).AND.(ylat(i,j).le.85))).OR.      &
419           (((xlong(i,j).ge.286).AND.(xlong(i,j).lt.345)).AND. &
420           ((ylat(i,j).ge.75).AND.(ylat(i,j).le.79.5))).OR.      &
421           (((xlong(i,j).ge.300).AND.(xlong(i,j).le.345)).AND. &
422           ((ylat(i,j).ge.67).AND.(ylat(i,j).le.75))).OR.    &
423           (((xlong(i,j).ge.305).AND.(xlong(i,j).le.330)).AND. &
424           ((ylat(i,j).ge.55).AND.(ylat(i,j).le.67))) ) mask_cal(i,j,2)=.true.
425   enddo
426enddo
427
428
429end subroutine input_clim
430!--------------------------------------------------------------------------------
431!subroutine input_climat_ref
432! quand on traite en absolu, pas besoin du climat de reference
433
434!end subroutine input_climat_ref
435
436
437
438SUBROUTINE init_forclim
439
440! fichiers snapshots
441  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
442
443! forcage temporel
444!------------------
445  namelist/forc_temporel/file_temporel,mincoefbmelt,maxcoefbmelt
446
447! lecture par namelist
448!---------------------
449! formats pour les ecritures dans 42
450428 format(A)
451  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
452  read(num_param,snap_forcage_mois_JB)
453  write(num_rep_42,428)'!___________________________________________________________' 
454  write(num_rep_42,428) '&snap_forcage_mois_JB                            ! module climat-forcage-insolation_mod'
455  write(num_rep_42,'(A,A)')   'filtr_t = ', filtr_t
456  write(num_rep_42,'(A,A)')   'filtr_p = ', filtr_p
457  write(num_rep_42,'(A,2(f7.1,","))')   'Seuil_haut = ', Seuil_haut(:)
458  write(num_rep_42,'(A,2(f7.1,","))')   'Seuil_bas = ', Seuil_bas(:)
459  write(num_rep_42,'(A,2(f5.1,","))')   'summorb = ', summorb(:)
460  write(num_rep_42,'(A,2(f7.1,","))')   'palier_ice = ', palier_ice(:,:)
461  write(num_rep_42,'(A,A)')   'surf_ice = ', surf_ice
462  write(num_rep_42,'(A,2(f3.1,","))')   'palier_CO2 = ', palier_CO2(:)
463  write(num_rep_42,'(A,A)')   'orb_file = ', orb_file
464  write(num_rep_42,'(A,A)')   'co2_file = ', co2_file
465  write(num_rep_42,*)'/'                     
466  write(num_rep_42,428) '! fichiers temperature et precip : 12 mois ordre : '
467  write(num_rep_42,428) '! 1/ CO2, 2/ regrowth, 3/ insol. insolmax_regrowthmaw_co2max --> insolmin_noice_co2min '
468  write(num_rep_42,*)
469
470!  do i=1,ntr
471! glaciaire
472!     write(filtr_t(i),'(A,i3,A)') filtr_t1(1:32),int(ttr(i)),filtr_t1(36:50)
473!     write(filtr_p(i),'(A,i3,A)') filtr_p1(1:34),int(ttr(i)),filtr_p1(38:52)
474!     write(filtr_t(i),'(A)') filtr_t
475!     write(filtr_p(i),'(A)') filtr_p
476!  enddo
477
478
479call lect_lapserate_months                               ! lit les lasperate mensuels
480                                                         ! pour une version spatialisee ecrire une autre routine
481! fichiers donnant l'evolution temporelle
482! ---------------------------- ------------
483
484rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
485read(num_param,forc_temporel)
486
487write(num_rep_42,428)'!___________________________________________________________' 
488write(num_rep_42,428) '&forc_temporel                            ! module climat_forcage_mois_mod'
489write(num_rep_42,'(A,A)')   'file_temporel =', file_temporel
490write(num_rep_42,*)     'mincoefbmelt      =', mincoefbmelt 
491write(num_rep_42,*)     'maxcoefbmelt      =', maxcoefbmelt 
492write(num_rep_42,*)'/'                     
493write(num_rep_42,428) '!fichier forcage temporel pour snapshot'
494write(num_rep_42,*)
495
496
497end subroutine init_forclim
498!---------------------------------------------------------------------
499
500!forcage climatique au cours du temps
501
502subroutine forclim
503
504implicit none
505
506real COEFT,COEFP       !
507!integer l              ! dumm index for loops on snapshots files  l=ITR,NTR-1
508!cdc integer itr            ! index of the current snapshot file (change with time)
509
510integer mo
511integer :: ictr
512integer :: igtr
513
514
515!Définition des variables propres au calcul de la fonction d'Insolation
516!y2_ins est le tableau contenant l'ensemble des dérivées secondes au points considérés (cf calcul de
517!l'interpolation par spline cubique).
518double precision,dimension(nb_points_laskar) :: y2_ins
519!insolval est la valeur de l'insolation d'été au pas de temps considere
520double precision :: insolval
521
522!initialisation des paramètres pour l'interpolation entre les fichiers de glace
523!real,allocatable :: sub_H(:,:)
524!real,allocatable :: sub_maillage_prox(:,:)
525!real :: nb_points,prox
526real :: prox
527!real,dimension(nx,ny) :: coeff1,coeff2
528real,dimension(nx,ny) :: coeff2
529integer :: k1,k2,m1,m2
530integer,dimension(nx,ny) :: Mask,Mask2
531integer :: update_H
532integer :: i,j,l
533
534!Definition des variables pour le calcul de la fonction de CO2
535double precision :: coefCO2
536
537!*****************************************************************************
538! triple interpolation orbite, glace, CO2 + correction altitude
539!*****************************************************************************
540
541
542select case(trim(CO2_function))
543case('constant')
544   CO2SnapMin=1
545   CO2SnapMax=1
546
547case('variable')
548   if (time.le.CO2_time(1)) then
549      CO2_value=CO2_val(1)
550   else if (time.ge.CO2_time(nb_slope_dis_CO2)) then
551      CO2_value=CO2_val(nb_slope_dis_CO2)
552   else
553      l=1
554      do while (time.gt.CO2_time(l))
555         l = l+1
556      end do
557     
558      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))
559 
560! pour ne pas depasser les valeurs CO2 des 2 snapshots min et max CO2
561      if (CO2_value.lt.palier_CO2(ctr)) then
562         CO2_value=palier_CO2(ctr)
563      else if (CO2_value.gt.palier_CO2(1)) then
564         CO2_value=palier_CO2(1)
565      endif
566
567   endif
568
569   l=1
570   do while (CO2_value.lt.palier_CO2(l+1))
571      l=l+1
572   enddo
573   CO2SnapMin=l
574   CO2SnapMax=l+1
575   print *,'time = ',time,'CO2_value = ',CO2_value, 'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax
576
577!   if (CO2_value.ge.360) then
578!      CO2SnapMin=1
579!      CO2SnapMax=2
580!   else if ((coefCO2.lt.360).and.(coefCO2.ge.280)) then
581!      CO2SnapMin=2
582!      CO2SnapMax=3
583!   else if (coefCO2.lt.280) then
584!      CO2SnapMin=3
585!      CO2SnapMax=4
586!   endif
587
588   !print *,'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax
589end select
590
591!*************************************************************************
592!*     triple interpolation orbite, glace, CO2 + correction altitude     *
593!*************************************************************************
594
595
596!Nouvelle fonction d'interpolation des param orbitaux qui permet de récupérer
597!l'insolation réelle de Laskar 2004
598
599
600select case(trim(orbital_interpolation))
601case('La2004')
602   
603   call spline(TEMPS_mod,INSOL,nb_points_laskar,yp1_ins,ypn_ins,y2_ins)
604
605   call splint(TEMPS_mod,INSOL,y2_ins,nb_points_laskar,insolval)
606 
607! Dans massud_param_list.dat, les snapshots doivent etre dans l'ordre
608! decroissant partant du plus chaud
609   if (insolval.gt.summorb(1)) then
610      Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,1)
611      Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,1)
612   else if (insolval.lt.summorb(ntr)) then
613      Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,ntr)
614      Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,ntr)
615   else
616      l=1
617      do while (insolval.lt.summorb(l))
618         l = l+1
619      enddo
620      factsin=(insolval-summorb(l))/(summorb(l-1)-summorb(l))
621      Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,l-1)*factsin+Tm_fin(:,:,:,:,:,l)*(1-factsin)
622      Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,l-1)*factsin+Pm_fin(:,:,:,:,:,l)*(1-factsin)
623   endif
624
625
626case('constant')
627   !Tm_time_fin_1(:,:,:,:,:)=(Tm_fin(:,:,:,:,:,1)+Tm_fin(:,:,:,:,:,2))/2
628   !Pm_time_fin_1(:,:,:,:,:)=(Pm_fin(:,:,:,:,:,1)+Pm_fin(:,:,:,:,:,2))/2
629   Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,1)
630   Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,1)
631
632
633case('sinusoid')
634   orbite=(1+sin(2*pi*time/40000.))/2.
635   Tm_time_fin_1(:,:,:,:,:)=orbite*Tm_fin(:,:,:,:,:,1)+(1-orbite)*Tm_fin(:,:,:,:,:,2)
636   Pm_time_fin_1(:,:,:,:,:)=orbite*Pm_fin(:,:,:,:,:,1)+(1-orbite)*Pm_fin(:,:,:,:,:,2)
637
638end select
639
640
641!do ictr=CO2SnapMin,CO2SnapMax,1
642!   do igtr=1,gtr
643!      print *,surf_ice(ictr,igtr)
644!   enddo
645!enddo
646
647
648!Correction d'altitude pour les états
649!il faut que Ssnap soit un tableau avec les différentes hauteurs de glace selon la
650!simulation (no ice, med ice, full ice et 2x, 2.5x, 3x)
651do j=1,ny
652   do i=1,nx
653      Zs(i,j)=max(slv(i,j),S(i,j))
654      do ictr=CO2SnapMin,CO2SnapMax,1
655         do igtr=1,gtr
656            if (Zs(i,j).ge.Ssnap(i,j,ictr,igtr)) then
657               do mo=1,mois
658                  ! correction d'altitude
659                  Tm_surf_mod(i,j,mo,ictr,igtr)=-lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,ictr,igtr)) &
660                                            +Tm_time_fin_1(i,j,mo,ictr,igtr)
661
662                  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) &
663                                            -Tm_time_fin_1(i,j,mo,ictr,igtr)))
664
665               end do
666            else if (Zs(i,j).lt.Ssnap(i,j,ictr,igtr)) then 
667               do mo=1,mois
668                  ! correction d'altitude avec condition pour que T ne soit
669                  ! pas supérieure à la T de l'état sans glace corrigé
670                  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), &
671                                            -lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,ictr,gtr))+Tm_time_fin_1(i,j,mo,ictr,gtr))
672
673                  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) &
674                                            -Tm_time_fin_1(i,j,mo,ictr,igtr)))
675               end do
676            endif
677         end do
678      end do
679   end do
680end do
681!print*,'debug Tm_surf_mod corr vert ', Tm_surf_mod(25,198,1,:,:)
682
683! calcul du nbr de points glace pose :
684do l=1,2
685   INP(l) = count(mask_cal(:,:,l).and.(H(:,:).gt.2.).and..not.flot(:,:))
686enddo
687np = INP(2)
688! np = count((H(:,:).gt.2.).and..not.flot(:,:)) ! version std
689!Interpolation glace
690surf_test=max(0,np)
691Height_new=H
692print *,'time = ',time
693!print *,'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax
694!print *,'Seuil_haut(CO2SnapMin) = ',Seuil_haut(CO2SnapMin),'Seuil_haut(CO2SnapMax) = ',Seuil_haut(CO2SnapMax)
695!print *,'Seuil_bas(CO2SnapMin) = ',Seuil_bas(CO2SnapMin),'Seuil_bas(CO2SnapMax) = ',Seuil_bas(CO2SnapMax)
696print *,'surf_test = ',surf_test
697!print *,'coeff1_stored(154,154) = ',coeff1_stored(154,154),'coeff2_stored(154,154) = ',coeff2_stored(154,154)
698!print *,'coeff1_stored(184,184) = ',coeff1_stored(184,184),'coeff2_stored(184,184) = ',coeff2_stored(184,184)
699!print *,'coeff1_stored(200,117) = ',coeff1_stored(200,117),'coeff2_stored(200,117) = ',coeff2_stored(200,117)
700
701select case(trim(ice_int))
702case('no')
703   Tm_time_fin_2(:,:,:,:)=Tm_surf_mod(:,:,:,:,1)
704   Pm_time_fin_2(:,:,:,:)=Pm_surf_mod(:,:,:,:,1)
705
706case('yes')
707   do ictr=CO2SnapMin,CO2SnapMax,1
708      if (surf_test.ge.Seuil_haut(ictr)) then
709         Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,1)
710         Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,1)
711
712         update_H=1
713      else if (surf_test.le.Seuil_bas(ictr)) then
714         Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,gtr)
715         Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,gtr)
716         
717         update_H=0
718      else
719         l=1
720         do while (surf_test.lt.palier_ice(ictr,l))
721            l = l+1
722         end do
723     
724         select case(trim(interpolation_ice))
725         case('linear')
726            retroicecoeff=(surf_test-palier_ice(ictr,l))/((palier_ice(ictr,l-1)-palier_ice(ictr,l)))
727            Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,l-1)*retroicecoeff+Tm_surf_mod(:,:,:,ictr,l)*(1-retroicecoeff)
728            Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,l-1)*retroicecoeff+Pm_surf_mod(:,:,:,ictr,l)*(1-retroicecoeff)
729
730         case('nonlinear')
731            !Important: the distinction between ice and vegetation points used
732            !before has been removed because it would end up with giving the
733            !same interpolation. In this code, the calculation is done only over
734            !surrounding ice points.
735            where (((Height_new > 5.).and.(Height_old < 5.)).or.((Height_new < 5.).and.(Height_old > 5.)))
736               Mask = 1.
737            elsewhere
738               Mask = 0.
739            end where
740           
741            !print *,'Mask = ',Mask
742
743!            allocate(sub_maillage_prox(2*indice+1,2*indice+1),stat=error)
744!            if (error.ne.0) then
745!               print *,'error: could not allocate memory for array'
746!            endif
747
748            Mask2 = 0.
749            do j=1,ny
750               do i=1,nx
751                  if (Mask(i,j).eq.1) then
752                     !print *,'i = ',i,'j = ',j,'Mask = ',Mask(i,j)
753                     do m2=max(1,-grid_int),min(grid_int,ny)  !cdc
754                        do m1=max(-grid_int,1),min(grid_int,nx)  !cdc
755                           Mask2(i+m1,j+m2) = 1.
756                        end do
757                     end do
758                  endif
759               end do
760            end do
761           
762            !print *,'Mask2 = ',Mask2
763            Hlarge(1:nx,1:ny)=H(:,:)
764            do j=1,ny
765               do i=1,nx
766                  if (Mask2(i,j).eq.1) then
767                     !The number of ice points nb_points is initialized at 0
768                     !Same for the proximity parameter, prox.
769                     !nb_points=0
770                     prox=0
771                     !Case where indice = 0, on the edges of the grid.
772                     !We actually dont care here because Antarctica is far from the
773                     !edges of the grid.
774                 
775                     !if ((j.eq.154).and.(i.eq.154)) then
776                     !   print *,'indice =',indice(i,j)
777                     !endif
778
779                     !If indice = 0, no interpolation, the climate used comes from
780                     !the deglaciated simulation (practically, we dont care about
781                     !this because it concerns only the 1st and last rows and
782                     !columns of the grid, so not Antarctica.... Could be
783                     !problematic in other contexts)
784                     !if (indice(i,j).eq.0) then
785                     !   Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)
786                     !   Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)
787                       
788                     !General case
789                     !else   
790                        !Allocation of the array sub_H, which will be used to count
791                        !the number of ice points around the focused point.
792                        !allocate(sub_H(2*indice(i,j)+1,2*indice(i,j)+1),stat=error)
793                        !if (error.ne.0) then
794                        !   print *,'error: could not allocate memory for array'
795                        !endif
796                        !sub_H is initialized at 1 everywhere except on the
797                        !focused point which is not taken into account.
798                        !sub_H=1
799                        !sub_H(indice(i,j)+1,indice(i,j)+1)=0
800                     
801                        !where (maillage_prox(indice(i,j),:,:) == 0.)
802                        !   sub_H(:,:)=0.
803                        !end where
804
805                        !Calculation of the number of ice points present around the
806                        !focused one.
807                        !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)
808                     
809                        !if ((j.eq.117).and.(i.eq.200)) then
810                        !   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)
811                        !   print *,'time =',time,'nb_points =',nb_points
812                        !endif
813                     
814                        !deallocate(sub_H)
815
816                        !Allocation of the array sub_maillage_prox, which is
817                        !extracted from the array maillage_prox (calculated above
818                        !in input_clim) accordingly to the position of the focused
819                        !point (ie the value of indice).
820!                        allocate(sub_maillage_prox(2*indice(i,j)+1,2*indice(i,j)+1),stat=error)
821!                        if (error.ne.0) then
822!                           print *,'error: could not allocate memory for array'
823!                        endif
824                        !sub_maillage_prox is initialized at 0 everywhere.
825!                        sub_maillage_prox=0
826                     
827                        !sub_maillage_prox is then filled with the corresponding
828                        !values from maillage_prox.
829!                        do k2=1,2*indice+1
830!                           do k1=1,2*indice+1
831!                              sub_maillage_prox(k1,k2)=maillage_prox(k1,k2)
832!                           end do
833!                        end do
834                        prox=0.
835                        !Calculation of the proximity parameter.
836!                        do k2=max(j-indice,1),min(j+indice,ny)
837!                           do k1=max(i-indice,1),min(i+indice,nx)
838!                              if (H(k1,k2).ge.5) prox=prox+sub_maillage_prox(k1,k2)
839!                           enddo
840!                        enddo
841                        prox=sum(maillage_prox,mask=Hlarge(i-grid_int:i+grid_int:1,j-grid_int:j+grid_int:1).ge.5)
842                       
843!                        if ((j.eq.198).and.(i.eq.25)) then
844!                           print *,'maillage_prox =',maillage_prox(:,:),'prox =',prox
845!                           print *,'prox_max =',prox_max(i,j)
846!                        endif
847   
848!                        deallocate(sub_maillage_prox)
849                     
850                        !Calculation of the relative coefficient (<1) for the
851                        !number of points and the proximity.
852                        !coeff1(i,j)=nb_points/nb_points_max(indice(i,j))
853                        coeff2(i,j)=prox/prox_max(i,j)
854                        !if (((j.eq.154).and.(i.eq.154)).or.((j.eq.184).and.(i.eq.184)).or.((j.eq.117).and.(i.eq.200))) then
855                        !       print *,'time = ',time,'i = ',i,'j = ',j,'coeff1 = ',coeff1(i,j),'coeff2 = ',coeff2(i,j)
856                        !endif
857
858                        !Interpolation:
859                        !Tm_surf_mod_noice*((1-a*coeff1)+(1-b*coeff2))/2+Tm_surf_mod_ice*(a*coeff1+b*coeff2)/2
860                        !Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)*((1-(2*poids_nb_points)*coeff1(i,j))+(1-(2*poids_proximity)*coeff2(i,j)))/2+Tm_surf_mod(i,j,:,ictr,l-1)*((2*poids_nb_points)*coeff1(i,j)+(2*poids_proximity)*coeff2(i,j))/2
861                        !Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)*((1-(2*poids_nb_points)*coeff1(i,j))+(1-(2*poids_proximity)*coeff2(i,j)))/2+Pm_surf_mod(i,j,:,ictr,l-1)*((2*poids_nb_points)*coeff1(i,j)+(2*poids_proximity)*coeff2(i,j))/2
862                       
863                        !Pour après
864                        Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)*(1-coeff2(i,j))+Tm_surf_mod(i,j,:,ictr,l-1)*coeff2(i,j)
865                        Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)*(1-coeff2(i,j))+Pm_surf_mod(i,j,:,ictr,l-1)*coeff2(i,j)
866                       
867                        !coeff1_stored(i,j)=coeff1(i,j)
868                        coeff2_stored(i,j)=coeff2(i,j)
869                        !print *,'coeff1_stored(154,154) = ',coeff1_stored(154,154),'coeff2_stored(154,154) = ',coeff2_stored(154,154)
870                        !print *,'coeff1_stored(184,184) = ',coeff1_stored(184,184),'coeff2_stored(184,184) = ',coeff2_stored(184,184)
871                        !print *,'coeff1_stored(200,117) = ',coeff1_stored(200,117),'coeff2_stored(200,117) = ',coeff2_stored(200,117)
872                     !endif
873                  else
874                     !Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)*((1-(2*poids_nb_points)*coeff1_stored(i,j))+(1-(2*poids_proximity)*coeff2_stored(i,j)))/2+Tm_surf_mod(i,j,:,ictr,l-1)*((2*poids_nb_points)*coeff1_stored(i,j)+(2*poids_proximity)*coeff2_stored(i,j))/2
875                     !Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)*((1-(2*poids_nb_points)*coeff1_stored(i,j))+(1-(2*poids_proximity)*coeff2_stored(i,j)))/2+Pm_surf_mod(i,j,:,ictr,l-1)*((2*poids_nb_points)*coeff1_stored(i,j)+(2*poids_proximity)*coeff2_stored(i,j))/2
876                     !Pour après
877                     Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l)*(1-coeff2_stored(i,j))+Tm_surf_mod(i,j,:,ictr,l-1)*coeff2_stored(i,j)
878                     Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l)*(1-coeff2_stored(i,j))+Pm_surf_mod(i,j,:,ictr,l-1)*coeff2_stored(i,j)
879                  endif
880
881               end do
882            end do
883!            deallocate(sub_maillage_prox)
884         end select
885
886         update_H=1
887      endif
888   end do
889end select
890!print*,'debug Tm_time_fin_2 ', Tm_time_fin_2(25,198,1,:)
891!print*,'debug coeff2_stored ', coeff2_stored(25,198)
892!print*,'debug grid_int ', grid_int
893
894
895!do j=1,ny
896!   do i=1,nx
897!      if (coeff1_stored(i,j).gt.0) then
898!         print *,'i = ',i,'j = ',j,'coeff1_stored = ',coeff1_stored(i,j)
899!      endif
900!      if (coeff2_stored(i,j).gt.0) then
901!         print *,'i = ',i,'j = ',j,'coeff2_stored = ',coeff2_stored(i,j)
902!      endif
903!   end do
904!end do
905
906if (update_H.eq.0) then
907   Height_old=0.
908   print *,'No update of H'
909else if (update_H.eq.1) then
910   Height_old=Height_new
911   print *,'H updated'
912else
913   print *,'Ice sheet topography error. Exit'
914   stop
915endif
916!print *,'time = ',time,'TEMP(200,117) = ',Tm_time_fin_2(200,117,:,ictr),'PRECIP(200,117) = ',Pm_time_fin_2(200,117,:,ictr)
917
918!fonction de décroissance du CO2
919!CO2 max en 1er puis décroissance dans l'ordre
920
921select case(trim(CO2_function))
922case('constant')
923   Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,1)
924   Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,1)
925
926
927case DEFAULT
928   print *,CO2_value
929   if (CO2_value.ge.palier_CO2(1)) then
930      Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,1)
931      Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,1)
932   else if (CO2_value.le.palier_CO2(ctr)) then
933      Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,ctr)
934      Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,ctr)
935   else
936      l=1
937      do while (CO2_value.lt.palier_CO2(l))
938         l = l+1
939      end do
940      select case(trim(interpolation_CO2))
941      case('linear')
942         factCO2=(CO2_value-palier_CO2(l))/((palier_CO2(l-1)-palier_CO2(l)))
943      case('logarithmic')
944         factCO2=log(CO2_value/palier_CO2(l))/log(palier_CO2(l-1)/palier_CO2(l))
945      end select
946      Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,l-1)*factCO2+Tm_time_fin_2(:,:,:,l)*(1-factCO2)
947      Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,l-1)*factCO2+Pm_time_fin_2(:,:,:,l)*(1-factCO2)
948   endif
949end select
950!print*,'debug Tm_time_fin_3 ', Tm_time_fin_3(25,198,1)
951!print *,time,h(196,145),uxbar(196,145),uybar(196,145),bm(196,145)
952
953do i=1,nx
954   do j=1,ny
955      do mo=1,mois
956         Tm_surf(i,j,mo)=Tm_time_fin_3(i,j,mo)
957         Pm_surf(i,j,mo)=Pm_time_fin_3(i,j,mo)
958      end do
959   end do
960end do
961!print*,'debug Tm_surf = Tm_3 ', Tm_surf(25,198,1)
962
963!coefbmshelf coefficient pour la fusion basale sous les ice shelves
964coefbmshelf=1.0
965coefbmshelf=max(coefbmshelf,mincoefbmelt)
966coefbmshelf=min(coefbmshelf,maxcoefbmelt)
967
968
969
970
971!*********************************
972!Correction d'altitude à vérifier
973!*********************************
974
975
976!  Correction d'altitude pour la temperature y compris sur les lacs
977! Zs est l'altitude de la surface qu'elle soit mer, glace ou lac
978
979!do j=1,ny
980!   do i=1,nx
981!      ZS(I,J)=max(slv(i,j),S(I,J))
982!      do mo=1,mois
983!         Tm_surf(i,j,mo)= - lapserate(i,j,mo) * (Zs(i,j)-Ssnap(i,j)) & ! correction d'altitude T
984!              + Tm_time(i,j,mo)
985!
986!         Pm_surf(i,j,mo)= Pm_time(i,j,mo)*exp(0.05*(Tm_surf(i,j,mo)-Tm_time(i,j,mo)))
987!
988!      end do
989!   end do
990!end do
991
992
993do mo=1,mois
994   Tmois(:,:,mo)=Tm_surf(:,:,mo)
995enddo
996!print*,'debug Tmois ', Tmois(25,198,1)
997
998! calcul de Tann et Tjuly pour les sorties :
999Tann(:,:)=sum(Tmois,dim=3)/12.     ! moy annuelle
1000Tjuly(:,:)=Tmois(:,:,7)            ! temp juillet
1001!print*,'debug Tann ', Tann(25,198)
1002
1003acc(:,:)=sum(Pm_surf,dim=3,mask=Tmois < psolid) ! /12.
1004acc(:,:)=acc(:,:)*1000./ro
1005
1006! pas beau ATTENTION
1007! pb plantage zone spitzberg bord de grille !!!!!
1008!do j=180,190
1009!   do i=123,126
1010!      acc(i,j)=0.
1011!      print*,'bsoc(i,j)',i,j,bsoc(i,j)
1012!   enddo
1013!enddo
1014
1015END subroutine forclim
1016
1017
1018
1019!************************************************************************
1020! Numerical Recipes
1021! interpolation spline cubique
1022
1023!Récupéré grace à Christophe. Modifié 17.04.13 par JB
1024SUBROUTINE splint(xa,ya,y2a,n,y)
1025  INTEGER,intent(in) :: n
1026  double precision,dimension(n),intent(in) :: xa,y2a,ya
1027  double precision,intent(out) :: y
1028! Calculates the cubic spline interpolation.
1029! Given 2 arrays of dimension n, xa and ya, and y2a, the second derivative
1030! of the function ya at any of the n points, it computes the interpolation for
1031! the array y of dimension nmax.
1032! for example, if n is 10001 (e.g. 1e6 years computed with Laskar algorithm and
1033! a sampling step of 100 years), then nmax would be 1e6 because GRISLI has a
1034! sampling step of 1 year.
1035  INTEGER l
1036  REAL a,b,h
1037! We will find the right place in the table by means of bisection.
1038
1039! do i=1,nmax
1040! print *,x(i)
1041! write (*,*) 'Press Enter to Continue'
1042! read (*,*)
1043! enddo
1044
1045!l=1
1046!do j=1,nmax
1047!if (j==1) then
1048!   l=2
1049!else if (j==nmax) then
1050!   l=n
1051!else if (x(j).gt.xa(l)) then
1052!   l=l+1
1053!endif
1054
1055
1056if (time==0) then
1057   l=2
1058else if (time==xa(n)) then
1059   l=n
1060else
1061   l=1
1062   do while (time.gt.xa(l))
1063   l=l+1
1064   enddo
1065endif
1066
1067
1068h=xa(l)-xa(l-1)
1069a=(xa(l)-time)/h
1070b=(time-xa(l-1))/h
1071y=a*ya(l-1)+b*ya(l)+((a**3-a)*y2a(l-1)+(b**3-b)*y2a(l))*(h**2)/6.
1072
1073! print *,time,l,xa(l-1),xa(l),ya(l-1),ya(l),h,a,b,y
1074! write (*,*) 'Press Enter to Continue'
1075! read (*,*)
1076
1077return
1078
1079END SUBROUTINE splint
1080
1081
1082
1083! Calculates the 2nd derivative of the y function for any points
1084! Recupere grace à Christophe. Modifié par JB 18.04.13
1085SUBROUTINE spline(x,y,n,yp1,ypn,y2)
1086  INTEGER,intent(in) :: n
1087  double precision,dimension(n), intent(in) :: x,y
1088  double precision,dimension(n), intent(out) :: y2
1089  double precision,intent(in) :: yp1,ypn
1090
1091  INTEGER i,k
1092  double precision :: p,qn,sig,un
1093  double precision,dimension(n) :: u
1094
1095if (yp1.gt..99e30) then ! The lower boundary condition is set either to 0
1096   y2(1)=0.
1097   u(1)=0.
1098else ! or else to have a specified first derivative.
1099   y2(1)=-0.5
1100   u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
1101endif
1102do i=2,n-1
1103! This is the decomposition loop of the tridiagonal algorithm. y2 and u are used
1104! for temporary storage of the decomposed factors.
1105
1106   sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
1107   p=sig*y2(i-1)+2.
1108   y2(i)=(sig-1.)/p
1109   u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
1110        /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
1111enddo
1112
1113if (ypn.gt..99e30) then ! The upper boundary condition is set eith
1114   qn=0.
1115   un=0.
1116else ! or else to have a specified first derivative.
1117   qn=0.5
1118   un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
1119endif
1120y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
1121do k=n-1,1,-1 ! This is the backsubstitution loop of the tridiagonal algorithm
1122   y2(k)=y2(k)*y2(k+1)+u(k)                             
1123enddo
1124return
1125
1126END SUBROUTINE spline
1127
1128
1129
1130
1131
1132
1133
1134!------------------------------------------------------------------------------------------------------------
1135subroutine lect_lapserate_months  ! lapserates mensuels mais uniformes spatialement
1136
1137implicit none
1138real,dimension(12) :: lect_lapse       ! pour la lecture
1139integer :: i,j
1140
1141namelist/lapse_month/lect_lapse
1142
1143! lecture de la namelist
1144
1145! formats pour les ecritures dans 42
1146428 format(A)
1147
1148rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
1149read(num_param,lapse_month)
1150
1151write(num_rep_42,428)'!___________________________________________________________' 
1152write(num_rep_42,428) '&lapse_month                                          !  module climat_forcage_mois_mod'
1153write(num_rep_42,'(A,12(f0.2,","))')  'lapse_month  = ', lect_lapse
1154write(num_rep_42,*)'/'                     
1155write(num_rep_42,428) '! laspe rates janvier -> decembre en deg/km'
1156
1157
1158! pour repasser en deg/m et copier dans lapserate
1159
1160do j=1,ny
1161   do i=1,nx
1162      lapserate(i,j,:)=lect_lapse(:)/1000.
1163   end do
1164end do
1165end subroutine lect_lapserate_months
1166!--------------------------------------------------------------------------------------------------------
1167
1168end module climat_forcage_insolation_mod
Note: See TracBrowser for help on using the repository browser.