!DERNIERE VERSION DU SCRIPT D'INTERPOLATION (May 2014) v2.1 !Interpolation NL/lineaire rayon 200/400 km + orbital forcing + CO2 decrease interpo lin ou log module climat_forcage_insolation_mod ! forcage avec champs mensuels ! ce qui a trait aux lacs proglaciaires passe dans un module separe. ! les tests sur geoplace sont enleves et sont remplacés par les ! lectures de fichers appropries. ! fonctionne avec un index en co2 et des snapshots a differents tx de co2 ! nouvelle version avec liste des variables utilisées par le module ! C. Dumas 07/2015 use module3D_phy,only:nx,ny,S,H,flot,Tann,Tjuly,Tmois,acc,coefbmshelf,sealevel_2d,ro,num_param,num_rep_42,dirnameinp,time,pi,dx,xlong,ylat,bsoc use netcdf use io_netcdf_grisli implicit none ! 1=decalaration variables !------------------------- integer :: nft ! NFT est le nombre de lignes a lire dans le fichier contenant le forcage climatique integer :: np ! nbr de points englaces poses integer,parameter :: mois=12 integer,parameter :: ntr=2 ! nb de snapshots selon les paramètres orbitaux integer,parameter :: gtr=4 ! nb de snapshots selon l'état de la calotte integer,parameter :: ctr=4 ! nb de snapshots selon le CO2 real,dimension(ctr) :: Seuil_haut real,dimension(ctr) :: Seuil_bas integer :: CO2SnapMin,CO2SnapMax character(len=12) :: orbital_interpolation='La2004' !Choose between La2004, constant and sinusoid character(len=12) :: ice_int='yes' !Choose between yes and no character(len=12) :: interpolation_ice='nonlinear' !Choose between linear and nonlinear integer :: rayon_interpolation_ice=60 !(In km) Choose a multiple of nx. !real :: poids_nb_points=0.5 !real :: poids_proximity=0.5 character(len=15) :: CO2_function='variable' !Choose between variable and constant integer :: nb_slope_dis_CO2 ! nbr de ligne fichier CO2 real :: CO2_value=980. character(len=12) :: interpolation_CO2='logarithmic' !Choose between linear and logarithmic !integer :: itr ! index of snapshot real :: coefsin ! coef pour param orbitaux double precision :: factsin real :: orbite real :: coefCO2,ccoefCO2 real :: factCO2 real :: surf_test real :: retroicecoeff real,dimension(ntr) :: summorb !Contient la moyenne DJF ou JJA de l'insolation pour ! chacun des paramètres orbitaux. On a WSO = 485 W.m-² et CSO = 382.2 W.m-² real,dimension(ctr,gtr) :: palier_ice real,dimension(ctr) :: palier_CO2 real,dimension(nx,ny,mois,ntr) :: Tm ! temperature mensuelle de chaque tranche real,dimension(nx,ny,mois,ntr) :: Pm ! precipitation mensuelle real,dimension(nx,ny,mois,ctr,gtr,ntr) :: Tm_fin ! tableau d'interpolation stylé real,dimension(nx,ny,mois,ctr,gtr,ntr) :: Pm_fin ! tableau d'interpolation stylé !interpolation sur param orbitaux !Définition des variables propres au calcul de la fonction d'Insolation !Nom du fichier contenant les points de temps et d'insolation récupérés du site !de Laskar. integer :: w !nb de points récupérés par le calcul de Laskar (e.g. 1e6 années avec un point !tous les 100 ans => 10001 points) integer,parameter :: nb_points_laskar=5001 !Les variables TEMPS, INSOL_dec, INSOL_jan, INSOL_fev correspondent au points !récupérés directement dans le fichier txt du site de Laskar. INSOL est la !moyenne des trois variables précédentes TEMPS_mod est une échelle de temps !recalibrée pour débuter à 0 et être échelonnée en années. double precision,dimension(nb_points_laskar) :: TEMPS,TEMPS_mod double precision,dimension(nb_points_laskar) :: INSOL,INSOL_dec,INSOL_jan,INSOL_fev !yp1_ins et ypn_ins sont les dérivés premières au 1er et dernier point de INSOL. double precision :: yp1_ins,ypn_ins !var_bidon1 et 2 sont des variables inutiles qui servent à ne pas multiplier le !nombre de tableaux contenant les différents points de TEMPS. real :: var_bidon1,var_bidon2 real,dimension(nx,ny,mois,ctr,gtr) :: Tm_time_fin_1 ! temperature mensuelle au temps time (non corrige altitude) real,dimension(nx,ny,mois,ctr,gtr) :: Pm_time_fin_1 ! precipitation mensuelle au temps time !interpolation sur surface glace !real,allocatable :: nb_points_max(:) real,dimension(nx,ny) :: prox_max ! somme des distances pour chaque point real,allocatable :: maillage_prox(:,:) ! distance des points autour du point etudie real,allocatable :: Hlarge(:,:) ! H avec grille elargie integer :: error real :: distsq real :: real_rayon_ice integer :: grid_int ! rayon de recherche en nbr de points de grille !integer :: indice real,dimension(nx,ny) :: Height_new,Height_old real,dimension(nx,ny) :: coeff2_stored !real,dimension(nx,ny) :: coeff1_stored,coeff2_stored real,dimension(nx,ny,mois,ctr) :: Tm_time_fin_2 ! temperature mensuelle au temps time (non corrige altitude) real,dimension(nx,ny,mois,ctr) :: Pm_time_fin_2 ! precipitation mensuelle au temps time !interpolation sur le CO2 double precision,allocatable :: CO2_time(:), CO2_val(:) ! depend du nbr de points du fichier CO2 real,dimension(nx,ny,mois) :: Tm_time_fin_3 ! temperature mensuelle au temps time (non corrige altitude) real,dimension(nx,ny,mois) :: Pm_time_fin_3 ! precipitation mensuelle au temps time real,dimension(nx,ny,mois,ctr,gtr) :: Tm_surf_mod ! correction altitude real,dimension(nx,ny,mois,ctr,gtr) :: Pm_surf_mod ! correction altitude real,dimension(nx,ny,mois) :: Tm_surf ! surface temperature (after topo. correction) real,dimension(nx,ny,mois) :: Pm_surf ! surface precipitation (after topo. correction) real,dimension(nx,ny,mois) :: lapserate ! lapse rate real :: psolid=2. ! temp limit between liquid and solid precip character(len=150) :: filin ! nom temporaire character(len=100) :: file_temporel ! forcage temporel character(len=100),dimension(ctr,gtr,ntr) :: filtr_t !snapshot temp file name character(len=100),dimension(ctr,gtr,ntr) :: filtr_p ! precip file character(len=100),dimension(ctr,gtr) :: surf_ice ! topo file character(len=100),dimension(3) :: orb_file ! orbite file character(len=100) :: co2_file ! co2 file real,dimension(nx,ny,ctr,gtr) :: Ssnap ! topo reference real,dimension(nx,ny) :: ZS !< surface topography above sea level ! specifique greeneem15 : calcul points englace sur Groenland uniquement : logical, dimension(nx,ny,2) :: mask_cal ! masque calotte Groenland (comme dans greeneem/output) integer,dimension(2) :: INP ! nbr de points englacés real :: mincoefbmelt ! butoirs pour coefbmshelf real :: maxcoefbmelt contains ! 2=lecture des inputs !-------------------- subroutine input_clim ! routine qui permet d'initialiser les variables climatiques ! variables locales !------------------- implicit none integer :: intr integer :: igtr integer :: ictr integer :: d1,d2 integer :: i,j ! character(len=100) :: file_ncdf !< fichier netcdf real*8, dimension(:,:,:), pointer :: data_3D => null() ! donnees lues dans le netcdf real*8, dimension(:,:),pointer :: data_2D => null() ! donnees lues dans le netcdf ! lecture des fichiers snapshots pour tout geoplace ! ------------------------------------------------- write(6,*) 'fichiers snapshots' do ictr=1,ctr ! co2 do igtr=1,gtr ! calotte do intr=1,ntr ! orbitaux !temperature filin=trim(dirnameinp)//'forcing/'//trim(filtr_t(ictr,igtr,intr)) call Read_ncdf_var('t2m',trim(filin),data_3D) ! Temperature Tm_fin(:,:,:,ictr,igtr,intr)=data_3D(:,:,:) ! write(6,*) trim(filin) ! open(20,file=trim(filin)) ! do j=1,ny ! do i=1,nx ! do J=ny,1,-1 ! do I=1,nx ! read(20,*) ti, tj, (Tm_fin(i,j,mo,ictr,igtr,intr),mo=1,12) ! end do ! end do ! close(20) !precipitation filin=trim(dirnameinp)//'forcing/'//trim(filtr_p(ictr,igtr,intr)) call Read_ncdf_var('precip',trim(filin),data_3D) ! precipitation Pm_fin(:,:,:,ictr,igtr,intr)=data_3D(:,:,:) ! write(6,*) trim(filin) ! open(20,file=trim(filin)) ! do j=1,ny ! do i=1,nx ! do J=ny,1,-1 ! do I=1,nx ! read(20,*) ti, tj, (Pm_fin(i,j,mo,ictr,igtr,intr),mo=1,12) ! end do ! end do ! close(20) end do ! topo filin=trim(dirnameinp)//'forcing/'//trim(surf_ice(ictr,igtr)) call Read_ncdf_var('TOPO',trim(filin),data_2D) ! topo Ssnap(:,:,ictr,igtr)=data_2D(:,:) where(Ssnap(:,:,ictr,igtr).eq.0.0) ! Pour PLIOMIP niv marin=25m Ssnap(:,:,ictr,igtr)=25.0 endwhere ! write(6,*) trim(filin) ! open(20,file=trim(filin)) ! do j=1,ny ! do i=1,nx ! do J=ny,1,-1 ! do I=1,nx ! read(20,*) Ssnap(i,j,ictr,igtr) ! end do ! end do ! close(20) end do end do !Lecture des parametres orbitaux select case(trim(orbital_interpolation)) case('La2004') print *,'Laskar 2004 orbital parameters' open(70,file=trim(dirnameinp)//'forcing/'//trim(orb_file(1))) do w=1,nb_points_laskar read(70,*) TEMPS(w),INSOL_dec(w) end do close(70) open(71,file=trim(dirnameinp)//'forcing/'//trim(orb_file(2))) do w=1,nb_points_laskar read(71,*) var_bidon1,INSOL_jan(w) end do close(71) open(72,file=trim(dirnameinp)//'forcing/'//trim(orb_file(3))) do w=1,nb_points_laskar read(72,*) var_bidon2,INSOL_fev(w) end do close(72) do w=1,nb_points_laskar TEMPS_mod(w) = 1000*TEMPS(w) !+34000000 INSOL(w)=(INSOL_dec(w)+INSOL_jan(w)+INSOL_fev(w))/3 end do !derivée 1er et dernier point yp1_ins=(INSOL(2)-INSOL(1))/(TEMPS_mod(2)-TEMPS_mod(1)) ypn_ins=(INSOL(nb_points_laskar)-INSOL(nb_points_laskar-1))/(TEMPS_mod(nb_points_laskar)-TEMPS_mod(nb_points_laskar-1)) case('constant') print *,'Constant orbital parameters' case('sinusoid') print *,'Sinusoidal variation of orbital parameters' end select !Calcul des poids relatifs lies a la position d'un point par rapport a un autre select case(trim(ice_int)) case('no') print *,"Ice interpolation not prescribed" case('yes') print *,"Ice interpolation prescribed" select case(trim(interpolation_ice)) case('linear') print *,"Linear ice mode" case('nonlinear') print *,"Non linear ice mode" ! Calculation of the radius of search in terms of indexes. grid_int=nint(rayon_interpolation_ice/(dx/1000.)) real_rayon_ice=real(rayon_interpolation_ice) print *,grid_int,ntr,gtr,ctr,intr,igtr,ictr print*,'dx=',dx print*,'rayon_interpolation_ice=',rayon_interpolation_ice print*,'grid_int=',grid_int ! indice=grid_int ! do j=1,ny ! do i=1,nx !Calculation of the array indice, which stores the number of points !taken into account in the interpolation. Practically overAntarctica, !it is always equal to the number of points within the circle defined !by grid_int, but closer to the edges of the grid, this number is smaller. !cdc indice(i,j)=min(min(i,j),min(nx-(i-1),ny-(j-1)))-1 !cdc if (indice(i,j).gt.(grid_int)) then ! indice(i,j)=grid_int !cdc endif ! if ((j.eq.198).and.(i.eq.25)) then ! print *,'time = ',time,'ictr = ',ictr,'np = ',np ! print *,'i = ',i,'j = ',j,'indice = ',indice(i,j) ! endif ! if ((j.eq.197).and.(i.eq.25)) then ! print *,'time = ',time,'ictr = ',ictr,'np = ',np ! print *,'i = ',i,'j = ',j,'indice = ',indice(i,j) ! endif ! enddo ! enddo !Allocation of the array nb_points_max, which stores the maximum number of !points potentially used in the interpolation for every indice value. !allocate(nb_points_max(grid_int),stat=error) !if (error.ne.0) then ! print *,'error: could not allocate memory for array' ! stop !endif !Allocation of the array prox_max. Same as above but for the proximity !parameter. !cdc allocate(prox_max(grid_int),stat=error) ! if (error.ne.0) then ! print *,'error: could not allocate memory for array' ! stop ! endif !Allocation of the array maillage_prox, which stores, for every indice !value, the proximity parameter associated with each point within the !radius of search. allocate(maillage_prox(-grid_int:grid_int,-grid_int:grid_int),stat=error) if (error.ne.0) then print *,'error: could not allocate memory for array' stop endif ! print*,'Hlarge',grid_int, 1-grid_int, nx+grid_int,ny+grid_int allocate(Hlarge(1-grid_int:nx+grid_int,1-grid_int:ny+grid_int),stat=error) if (error.ne.0) then print *,'error: could not allocate memory for array' stop endif Hlarge=0. !Initialization at 0 for the three arrays allocated before. !nb_points_max=0 prox_max=0 maillage_prox=0 !Filling of the three arrays. The proximity parameter depends on a ration !of exponential functions. !cdc do k=1,grid_int do d2=-grid_int,grid_int do d1=-grid_int,grid_int distsq=(abs(d1)*(dx/1000.))**2+(abs(d2)*(dx/1000.))**2 !print *,'distsq = ',distsq,'real_rayon_ice =',real_rayon_ice if ((sqrt(distsq).le.real_rayon_ice).and.((d1/=0).or.(d2/=0))) then !nb_points_max(k)=nb_points_max(k)+1 !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.)) !Indexes displacement to match elements 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.)) !print *,'nb_points_max = ',nb_points_max(k),'prox_max = ',prox_max(k),'maillage_prox = ',maillage_prox(k,d1+(k+1),d2+(k+1)) endif enddo enddo !cdc enddo prox_max(:,:)=SUM(maillage_prox(:,:)) ! PRINT *,'maillage_prox =',maillage_prox(:,:) ! PRINT*,'sum(maillage_prox)',SUM(maillage_prox(:,:)) ! PRINT*,'fin calcul maillage prox' end select end select Height_new=0. Height_old=0. !coeff1_stored=0. coeff2_stored=0. select case(trim(CO2_function)) case('constant') print *,'CO2 kept constant at ',CO2_value case('variable') print *,'The CO2 will vary with time' open(75,file=trim(dirnameinp)//'forcing/'//trim(co2_file)) read(75,*) nb_slope_dis_CO2 allocate(CO2_time(nb_slope_dis_CO2),stat=error) if (error.ne.0) then print *,'error: could not allocate memory for array CO2_time' stop endif allocate(CO2_val(nb_slope_dis_CO2),stat=error) if (error.ne.0) then print *,'error: could not allocate memory for array CO2_val' stop endif do w=1,nb_slope_dis_CO2 read(75,*) CO2_time(w),CO2_val(w) end do close(75) !print *,'CO2_time = ',CO2_time(:),'CO2_val = ',CO2_val(:) end select ! initilisation region Groenland pour calcul nbr points englaces mask_cal(:,:,:)=.false. mask_cal(:,:,1)=.true. do j=1,ny do i=1,nx ! -- Groenland IF ( (((xlong(i,j).ge.290).AND.(xlong(i,j).le.350)).AND. & ((ylat(i,j).ge.79.5).AND.(ylat(i,j).le.85))).OR. & (((xlong(i,j).ge.286).AND.(xlong(i,j).lt.345)).AND. & ((ylat(i,j).ge.75).AND.(ylat(i,j).le.79.5))).OR. & (((xlong(i,j).ge.300).AND.(xlong(i,j).le.345)).AND. & ((ylat(i,j).ge.67).AND.(ylat(i,j).le.75))).OR. & (((xlong(i,j).ge.305).AND.(xlong(i,j).le.330)).AND. & ((ylat(i,j).ge.55).AND.(ylat(i,j).le.67))) ) mask_cal(i,j,2)=.true. enddo enddo end subroutine input_clim !-------------------------------------------------------------------------------- !subroutine input_climat_ref ! quand on traite en absolu, pas besoin du climat de reference !end subroutine input_climat_ref SUBROUTINE init_forclim ! fichiers snapshots 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 ! forcage temporel !------------------ namelist/forc_temporel/file_temporel,mincoefbmelt,maxcoefbmelt ! lecture par namelist !--------------------- ! formats pour les ecritures dans 42 428 format(A) rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,snap_forcage_mois_JB) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&snap_forcage_mois_JB ! module climat-forcage-insolation_mod' write(num_rep_42,'(A,A)') 'filtr_t = ', filtr_t write(num_rep_42,'(A,A)') 'filtr_p = ', filtr_p write(num_rep_42,'(A,2(f7.1,","))') 'Seuil_haut = ', Seuil_haut(:) write(num_rep_42,'(A,2(f7.1,","))') 'Seuil_bas = ', Seuil_bas(:) write(num_rep_42,'(A,2(f5.1,","))') 'summorb = ', summorb(:) write(num_rep_42,'(A,2(f7.1,","))') 'palier_ice = ', palier_ice(:,:) write(num_rep_42,'(A,A)') 'surf_ice = ', surf_ice write(num_rep_42,'(A,2(f6.1,","))') 'palier_CO2 = ', palier_CO2(:) write(num_rep_42,'(A,A)') 'orb_file = ', orb_file write(num_rep_42,'(A,A)') 'co2_file = ', co2_file write(num_rep_42,*)'/' write(num_rep_42,428) '! fichiers temperature et precip : 12 mois ordre : ' write(num_rep_42,428) '! 1/ CO2, 2/ regrowth, 3/ insol. insolmax_regrowthmaw_co2max --> insolmin_noice_co2min ' write(num_rep_42,*) ! do i=1,ntr ! glaciaire ! write(filtr_t(i),'(A,i3,A)') filtr_t1(1:32),int(ttr(i)),filtr_t1(36:50) ! write(filtr_p(i),'(A,i3,A)') filtr_p1(1:34),int(ttr(i)),filtr_p1(38:52) ! write(filtr_t(i),'(A)') filtr_t ! write(filtr_p(i),'(A)') filtr_p ! enddo call lect_lapserate_months ! lit les lasperate mensuels ! pour une version spatialisee ecrire une autre routine ! fichiers donnant l'evolution temporelle ! ---------------------------- ------------ rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,forc_temporel) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&forc_temporel ! module climat_forcage_mois_mod' write(num_rep_42,'(A,A)') 'file_temporel =', file_temporel write(num_rep_42,*) 'mincoefbmelt =', mincoefbmelt write(num_rep_42,*) 'maxcoefbmelt =', maxcoefbmelt write(num_rep_42,*)'/' write(num_rep_42,428) '!fichier forcage temporel pour snapshot' write(num_rep_42,*) end subroutine init_forclim !--------------------------------------------------------------------- !forcage climatique au cours du temps subroutine forclim implicit none !integer l ! dumm index for loops on snapshots files l=ITR,NTR-1 !cdc integer itr ! index of the current snapshot file (change with time) integer mo integer :: ictr integer :: igtr !Définition des variables propres au calcul de la fonction d'Insolation !y2_ins est le tableau contenant l'ensemble des dérivées secondes au points considérés (cf calcul de !l'interpolation par spline cubique). double precision,dimension(nb_points_laskar) :: y2_ins !insolval est la valeur de l'insolation d'été au pas de temps considere double precision :: insolval !initialisation des paramètres pour l'interpolation entre les fichiers de glace !real,allocatable :: sub_H(:,:) !real,allocatable :: sub_maillage_prox(:,:) !real :: nb_points,prox real :: prox !real,dimension(nx,ny) :: coeff1,coeff2 real,dimension(nx,ny) :: coeff2 integer :: m1,m2 integer,dimension(nx,ny) :: Mask,Mask2 integer :: update_H integer :: i,j,l !Definition des variables pour le calcul de la fonction de CO2 !double precision :: coefCO2 !***************************************************************************** ! triple interpolation orbite, glace, CO2 + correction altitude !***************************************************************************** select case(trim(CO2_function)) case('constant') CO2SnapMin=1 CO2SnapMax=1 case('variable') if (time.le.CO2_time(1)) then CO2_value=CO2_val(1) else if (time.ge.CO2_time(nb_slope_dis_CO2)) then CO2_value=CO2_val(nb_slope_dis_CO2) else l=1 do while (time.gt.CO2_time(l)) l = l+1 end do 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)) ! pour ne pas depasser les valeurs CO2 des 2 snapshots min et max CO2 if (CO2_value.lt.palier_CO2(ctr)) then CO2_value=palier_CO2(ctr) else if (CO2_value.gt.palier_CO2(1)) then CO2_value=palier_CO2(1) endif endif l=1 do while (CO2_value.lt.palier_CO2(l+1)) l=l+1 enddo CO2SnapMin=l CO2SnapMax=l+1 print *,'time = ',time,'CO2_value = ',CO2_value, 'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax ! if (CO2_value.ge.360) then ! CO2SnapMin=1 ! CO2SnapMax=2 ! else if ((coefCO2.lt.360).and.(coefCO2.ge.280)) then ! CO2SnapMin=2 ! CO2SnapMax=3 ! else if (coefCO2.lt.280) then ! CO2SnapMin=3 ! CO2SnapMax=4 ! endif !print *,'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax end select !************************************************************************* !* triple interpolation orbite, glace, CO2 + correction altitude * !************************************************************************* !Nouvelle fonction d'interpolation des param orbitaux qui permet de récupérer !l'insolation réelle de Laskar 2004 select case(trim(orbital_interpolation)) case('La2004') call spline(TEMPS_mod,INSOL,nb_points_laskar,yp1_ins,ypn_ins,y2_ins) call splint(TEMPS_mod,INSOL,y2_ins,nb_points_laskar,insolval) ! Dans massud_param_list.dat, les snapshots doivent etre dans l'ordre ! decroissant partant du plus chaud if (insolval.gt.summorb(1)) then Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,1) Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,1) else if (insolval.lt.summorb(ntr)) then Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,ntr) Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,ntr) else l=1 do while (insolval.lt.summorb(l)) l = l+1 enddo factsin=(insolval-summorb(l))/(summorb(l-1)-summorb(l)) Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,l-1)*factsin+Tm_fin(:,:,:,:,:,l)*(1-factsin) Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,l-1)*factsin+Pm_fin(:,:,:,:,:,l)*(1-factsin) endif case('constant') !Tm_time_fin_1(:,:,:,:,:)=(Tm_fin(:,:,:,:,:,1)+Tm_fin(:,:,:,:,:,2))/2 !Pm_time_fin_1(:,:,:,:,:)=(Pm_fin(:,:,:,:,:,1)+Pm_fin(:,:,:,:,:,2))/2 Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,1) Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,1) case('sinusoid') orbite=(1+sin(2*pi*time/40000.))/2. Tm_time_fin_1(:,:,:,:,:)=orbite*Tm_fin(:,:,:,:,:,1)+(1-orbite)*Tm_fin(:,:,:,:,:,2) Pm_time_fin_1(:,:,:,:,:)=orbite*Pm_fin(:,:,:,:,:,1)+(1-orbite)*Pm_fin(:,:,:,:,:,2) end select !do ictr=CO2SnapMin,CO2SnapMax,1 ! do igtr=1,gtr ! print *,surf_ice(ictr,igtr) ! enddo !enddo !Correction d'altitude pour les états !il faut que Ssnap soit un tableau avec les différentes hauteurs de glace selon la !simulation (no ice, med ice, full ice et 2x, 2.5x, 3x) do j=1,ny do i=1,nx Zs(i,j)=max(sealevel_2d(i,j),S(i,j)) do ictr=CO2SnapMin,CO2SnapMax,1 do igtr=1,gtr if (Zs(i,j).ge.Ssnap(i,j,ictr,igtr)) then do mo=1,mois ! correction d'altitude Tm_surf_mod(i,j,mo,ictr,igtr)=-lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,ictr,igtr)) & +Tm_time_fin_1(i,j,mo,ictr,igtr) 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) & -Tm_time_fin_1(i,j,mo,ictr,igtr))) end do else if (Zs(i,j).lt.Ssnap(i,j,ictr,igtr)) then do mo=1,mois ! correction d'altitude avec condition pour que T ne soit ! pas supérieure à la T de l'état sans glace corrigé 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), & -lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,ictr,gtr))+Tm_time_fin_1(i,j,mo,ictr,gtr)) 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) & -Tm_time_fin_1(i,j,mo,ictr,igtr))) end do endif end do end do end do end do !print*,'debug Tm_surf_mod corr vert ', Tm_surf_mod(25,198,1,:,:) ! calcul du nbr de points glace pose : do l=1,2 INP(l) = count(mask_cal(:,:,l).and.(H(:,:).gt.2.).and..not.flot(:,:)) enddo np = INP(2) ! np = count((H(:,:).gt.2.).and..not.flot(:,:)) ! version std !Interpolation glace surf_test=max(0,np) Height_new=H print *,'time = ',time !print *,'CO2SnapMin = ',CO2SnapMin,'CO2SnapMax = ',CO2SnapMax !print *,'Seuil_haut(CO2SnapMin) = ',Seuil_haut(CO2SnapMin),'Seuil_haut(CO2SnapMax) = ',Seuil_haut(CO2SnapMax) !print *,'Seuil_bas(CO2SnapMin) = ',Seuil_bas(CO2SnapMin),'Seuil_bas(CO2SnapMax) = ',Seuil_bas(CO2SnapMax) print *,'surf_test = ',surf_test !print *,'coeff1_stored(154,154) = ',coeff1_stored(154,154),'coeff2_stored(154,154) = ',coeff2_stored(154,154) !print *,'coeff1_stored(184,184) = ',coeff1_stored(184,184),'coeff2_stored(184,184) = ',coeff2_stored(184,184) !print *,'coeff1_stored(200,117) = ',coeff1_stored(200,117),'coeff2_stored(200,117) = ',coeff2_stored(200,117) select case(trim(ice_int)) case('no') Tm_time_fin_2(:,:,:,:)=Tm_surf_mod(:,:,:,:,1) Pm_time_fin_2(:,:,:,:)=Pm_surf_mod(:,:,:,:,1) case('yes') do ictr=CO2SnapMin,CO2SnapMax,1 if (surf_test.ge.Seuil_haut(ictr)) then Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,1) Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,1) update_H=1 else if (surf_test.le.Seuil_bas(ictr)) then Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,gtr) Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,gtr) update_H=0 else l=1 do while (surf_test.lt.palier_ice(ictr,l)) l = l+1 end do select case(trim(interpolation_ice)) case('linear') retroicecoeff=(surf_test-palier_ice(ictr,l))/((palier_ice(ictr,l-1)-palier_ice(ictr,l))) Tm_time_fin_2(:,:,:,ictr)=Tm_surf_mod(:,:,:,ictr,l-1)*retroicecoeff+Tm_surf_mod(:,:,:,ictr,l)*(1-retroicecoeff) Pm_time_fin_2(:,:,:,ictr)=Pm_surf_mod(:,:,:,ictr,l-1)*retroicecoeff+Pm_surf_mod(:,:,:,ictr,l)*(1-retroicecoeff) case('nonlinear') !Important: the distinction between ice and vegetation points used !before has been removed because it would end up with giving the !same interpolation. In this code, the calculation is done only over !surrounding ice points. where (((Height_new > 5.).and.(Height_old < 5.)).or.((Height_new < 5.).and.(Height_old > 5.))) Mask = 1. elsewhere Mask = 0. end where !print *,'Mask = ',Mask ! allocate(sub_maillage_prox(2*indice+1,2*indice+1),stat=error) ! if (error.ne.0) then ! print *,'error: could not allocate memory for array' ! endif Mask2 = 0. do j=1,ny do i=1,nx if (Mask(i,j).eq.1) then !print *,'i = ',i,'j = ',j,'Mask = ',Mask(i,j) do m2=max(1,j-grid_int),min(j+grid_int,ny) !cdc do m1=max(i-grid_int,1),min(i+grid_int,nx) !cdc Mask2(m1,m2) = 1. end do end do endif end do end do !print *,'Mask2 = ',Mask2 Hlarge(1:nx,1:ny)=H(:,:) do j=1,ny do i=1,nx if (Mask2(i,j).eq.1) then !The number of ice points nb_points is initialized at 0 !Same for the proximity parameter, prox. !nb_points=0 ! prox=0 !Case where indice = 0, on the edges of the grid. !We actually dont care here because Antarctica is far from the !edges of the grid. !if ((j.eq.154).and.(i.eq.154)) then ! print *,'indice =',indice(i,j) !endif !If indice = 0, no interpolation, the climate used comes from !the deglaciated simulation (practically, we dont care about !this because it concerns only the 1st and last rows and !columns of the grid, so not Antarctica.... Could be !problematic in other contexts) !if (indice(i,j).eq.0) then ! Tm_time_fin_2(i,j,:,ictr)=Tm_surf_mod(i,j,:,ictr,l) ! Pm_time_fin_2(i,j,:,ictr)=Pm_surf_mod(i,j,:,ictr,l) !General case !else !Allocation of the array sub_H, which will be used to count !the number of ice points around the focused point. !allocate(sub_H(2*indice(i,j)+1,2*indice(i,j)+1),stat=error) !if (error.ne.0) then ! print *,'error: could not allocate memory for array' !endif !sub_H is initialized at 1 everywhere except on the !focused point which is not taken into account. !sub_H=1 !sub_H(indice(i,j)+1,indice(i,j)+1)=0 !where (maillage_prox(indice(i,j),:,:) == 0.) ! sub_H(:,:)=0. !end where !Calculation of the number of ice points present around the !focused one. !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) !if ((j.eq.117).and.(i.eq.200)) then ! 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) ! print *,'time =',time,'nb_points =',nb_points !endif !deallocate(sub_H) !Allocation of the array sub_maillage_prox, which is !extracted from the array maillage_prox (calculated above !in input_clim) accordingly to the position of the focused !point (ie the value of indice). ! allocate(sub_maillage_prox(2*indice(i,j)+1,2*indice(i,j)+1),stat=error) ! if (error.ne.0) then ! print *,'error: could not allocate memory for array' ! endif !sub_maillage_prox is initialized at 0 everywhere. ! sub_maillage_prox=0 !sub_maillage_prox is then filled with the corresponding !values from maillage_prox. ! do k2=1,2*indice+1 ! do k1=1,2*indice+1 ! sub_maillage_prox(k1,k2)=maillage_prox(k1,k2) ! end do ! end do prox=0. !Calculation of the proximity parameter. ! do k2=max(j-indice,1),min(j+indice,ny) ! do k1=max(i-indice,1),min(i+indice,nx) ! if (H(k1,k2).ge.5) prox=prox+sub_maillage_prox(k1,k2) ! enddo ! enddo ! maillage_prox_ij(i-grid_int:i+grid_int:1,j-grid_int:j+grid_int:1)=maillage_prox prox=sum(maillage_prox,mask=Hlarge(i-grid_int:i+grid_int:1,j-grid_int:j+grid_int:1).ge.5) ! if ((j.eq.198).and.(i.eq.25)) then ! print *,'maillage_prox =',maillage_prox(:,:),'prox =',prox ! print *,'prox_max =',prox_max(i,j) ! endif ! deallocate(sub_maillage_prox) !Calculation of the relative coefficient (<1) for the !number of points and the proximity. !coeff1(i,j)=nb_points/nb_points_max(indice(i,j)) coeff2(i,j)=prox/prox_max(i,j) ! PRINT*,'i=',i,'j=',j ! PRINT*,'prox=',prox ! PRINT*,'prox_max=',prox_max(i,j) ! PRINT*,'coeff2=',coeff2(i,j) ! DO m1=j+grid_int,j-grid_int,-1 ! PRINT*,'Hlarge',Hlarge(i-grid_int:i+grid_int:1,m1) ! ENDDO ! PRINT* !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 ! print *,'time = ',time,'i = ',i,'j = ',j,'coeff1 = ',coeff1(i,j),'coeff2 = ',coeff2(i,j) !endif !Interpolation: !Tm_surf_mod_noice*((1-a*coeff1)+(1-b*coeff2))/2+Tm_surf_mod_ice*(a*coeff1+b*coeff2)/2 !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 !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 !Pour après 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) 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) !coeff1_stored(i,j)=coeff1(i,j) coeff2_stored(i,j)=coeff2(i,j) !print *,'coeff1_stored(154,154) = ',coeff1_stored(154,154),'coeff2_stored(154,154) = ',coeff2_stored(154,154) !print *,'coeff1_stored(184,184) = ',coeff1_stored(184,184),'coeff2_stored(184,184) = ',coeff2_stored(184,184) !print *,'coeff1_stored(200,117) = ',coeff1_stored(200,117),'coeff2_stored(200,117) = ',coeff2_stored(200,117) !endif else !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 !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 !Pour après 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) 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) endif end do end do ! deallocate(sub_maillage_prox) end select update_H=1 endif end do end select !print*,'debug Tm_time_fin_2 ', Tm_time_fin_2(25,198,1,:) !print*,'debug coeff2_stored ', coeff2_stored(25,198) !print*,'debug grid_int ', grid_int !do j=1,ny ! do i=1,nx ! if (coeff1_stored(i,j).gt.0) then ! print *,'i = ',i,'j = ',j,'coeff1_stored = ',coeff1_stored(i,j) ! endif ! if (coeff2_stored(i,j).gt.0) then ! print *,'i = ',i,'j = ',j,'coeff2_stored = ',coeff2_stored(i,j) ! endif ! end do !end do if (update_H.eq.0) then Height_old=0. print *,'No update of H' else if (update_H.eq.1) then Height_old=Height_new print *,'H updated' else print *,'Ice sheet topography error. Exit' stop endif !print *,'time = ',time,'TEMP(200,117) = ',Tm_time_fin_2(200,117,:,ictr),'PRECIP(200,117) = ',Pm_time_fin_2(200,117,:,ictr) !fonction de décroissance du CO2 !CO2 max en 1er puis décroissance dans l'ordre select case(trim(CO2_function)) case('constant') Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,1) Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,1) case DEFAULT print *,CO2_value if (CO2_value.ge.palier_CO2(1)) then Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,1) Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,1) else if (CO2_value.le.palier_CO2(ctr)) then Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,ctr) Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,ctr) else l=1 do while (CO2_value.lt.palier_CO2(l)) l = l+1 end do select case(trim(interpolation_CO2)) case('linear') factCO2=(CO2_value-palier_CO2(l))/((palier_CO2(l-1)-palier_CO2(l))) case('logarithmic') factCO2=log(CO2_value/palier_CO2(l))/log(palier_CO2(l-1)/palier_CO2(l)) end select Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,l-1)*factCO2+Tm_time_fin_2(:,:,:,l)*(1-factCO2) Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,l-1)*factCO2+Pm_time_fin_2(:,:,:,l)*(1-factCO2) endif end select !print*,'debug Tm_time_fin_3 ', Tm_time_fin_3(25,198,1) !print *,time,h(196,145),uxbar(196,145),uybar(196,145),bm(196,145) do i=1,nx do j=1,ny do mo=1,mois Tm_surf(i,j,mo)=Tm_time_fin_3(i,j,mo) Pm_surf(i,j,mo)=Pm_time_fin_3(i,j,mo) end do end do end do !print*,'debug Tm_surf = Tm_3 ', Tm_surf(25,198,1) !coefbmshelf coefficient pour la fusion basale sous les ice shelves coefbmshelf=1.0 coefbmshelf=max(coefbmshelf,mincoefbmelt) coefbmshelf=min(coefbmshelf,maxcoefbmelt) !********************************* !Correction d'altitude à vérifier !********************************* ! Correction d'altitude pour la temperature y compris sur les lacs ! Zs est l'altitude de la surface qu'elle soit mer, glace ou lac !do j=1,ny ! do i=1,nx ! ZS(I,J)=max(slv(i,j),S(I,J)) ! do mo=1,mois ! Tm_surf(i,j,mo)= - lapserate(i,j,mo) * (Zs(i,j)-Ssnap(i,j)) & ! correction d'altitude T ! + Tm_time(i,j,mo) ! ! Pm_surf(i,j,mo)= Pm_time(i,j,mo)*exp(0.05*(Tm_surf(i,j,mo)-Tm_time(i,j,mo))) ! ! end do ! end do !end do do mo=1,mois Tmois(:,:,mo)=Tm_surf(:,:,mo) enddo !print*,'debug Tmois ', Tmois(25,198,1) ! calcul de Tann et Tjuly pour les sorties : Tann(:,:)=sum(Tmois,dim=3)/12. ! moy annuelle Tjuly(:,:)=Tmois(:,:,7) ! temp juillet !print*,'debug Tann ', Tann(25,198) acc(:,:)=sum(Pm_surf,dim=3,mask=Tmois < psolid) ! /12. acc(:,:)=acc(:,:)*1000./ro ! pas beau ATTENTION ! pb plantage zone spitzberg bord de grille !!!!! !do j=180,190 ! do i=123,126 ! acc(i,j)=0. ! print*,'bsoc(i,j)',i,j,bsoc(i,j) ! enddo !enddo END subroutine forclim !************************************************************************ ! Numerical Recipes ! interpolation spline cubique !Récupéré grace à Christophe. Modifié 17.04.13 par JB SUBROUTINE splint(xa,ya,y2a,n,y) INTEGER,intent(in) :: n double precision,dimension(n),intent(in) :: xa,y2a,ya double precision,intent(out) :: y ! Calculates the cubic spline interpolation. ! Given 2 arrays of dimension n, xa and ya, and y2a, the second derivative ! of the function ya at any of the n points, it computes the interpolation for ! the array y of dimension nmax. ! for example, if n is 10001 (e.g. 1e6 years computed with Laskar algorithm and ! a sampling step of 100 years), then nmax would be 1e6 because GRISLI has a ! sampling step of 1 year. INTEGER l REAL a,b,h ! We will find the right place in the table by means of bisection. ! do i=1,nmax ! print *,x(i) ! write (*,*) 'Press Enter to Continue' ! read (*,*) ! enddo !l=1 !do j=1,nmax !if (j==1) then ! l=2 !else if (j==nmax) then ! l=n !else if (x(j).gt.xa(l)) then ! l=l+1 !endif if (time==0) then l=2 else if (time==xa(n)) then l=n else l=1 do while (time.gt.xa(l)) l=l+1 enddo endif h=xa(l)-xa(l-1) a=(xa(l)-time)/h b=(time-xa(l-1))/h y=a*ya(l-1)+b*ya(l)+((a**3-a)*y2a(l-1)+(b**3-b)*y2a(l))*(h**2)/6. ! print *,time,l,xa(l-1),xa(l),ya(l-1),ya(l),h,a,b,y ! write (*,*) 'Press Enter to Continue' ! read (*,*) return END SUBROUTINE splint ! Calculates the 2nd derivative of the y function for any points ! Recupere grace à Christophe. Modifié par JB 18.04.13 SUBROUTINE spline(x,y,n,yp1,ypn,y2) INTEGER,intent(in) :: n double precision,dimension(n), intent(in) :: x,y double precision,dimension(n), intent(out) :: y2 double precision,intent(in) :: yp1,ypn INTEGER i,k double precision :: p,qn,sig,un double precision,dimension(n) :: u if (yp1.gt..99e30) then ! The lower boundary condition is set either to 0 y2(1)=0. u(1)=0. else ! or else to have a specified first derivative. y2(1)=-0.5 u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) endif do i=2,n-1 ! This is the decomposition loop of the tridiagonal algorithm. y2 and u are used ! for temporary storage of the decomposed factors. sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) p=sig*y2(i-1)+2. y2(i)=(sig-1.)/p u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) & /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p enddo if (ypn.gt..99e30) then ! The upper boundary condition is set eith qn=0. un=0. else ! or else to have a specified first derivative. qn=0.5 un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) endif y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) do k=n-1,1,-1 ! This is the backsubstitution loop of the tridiagonal algorithm y2(k)=y2(k)*y2(k+1)+u(k) enddo return END SUBROUTINE spline !------------------------------------------------------------------------------------------------------------ subroutine lect_lapserate_months ! lapserates mensuels mais uniformes spatialement implicit none real,dimension(12) :: lect_lapse ! pour la lecture integer :: i,j namelist/lapse_month/lect_lapse ! lecture de la namelist ! formats pour les ecritures dans 42 428 format(A) rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,lapse_month) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&lapse_month ! module climat_forcage_mois_mod' write(num_rep_42,'(A,12(f0.2,","))') 'lapse_month = ', lect_lapse write(num_rep_42,*)'/' write(num_rep_42,428) '! laspe rates janvier -> decembre en deg/km' ! pour repasser en deg/m et copier dans lapserate do j=1,ny do i=1,nx lapserate(i,j,:)=lect_lapse(:)/1000. end do end do end subroutine lect_lapserate_months !-------------------------------------------------------------------------------------------------------- end module climat_forcage_insolation_mod