!Interpolation NL/lineaire rayon 200/400 km + orbital forcing + CO2 decrease interpo lin ou !log ! Ajout de l'interpolation verticale sous maille LMDZ -> GRISLI C. DUMAS Fev 2015 ! Lecture d'un fichier de topo associé au fichier climat module climat_forcage_insolation_mod_oneway ! forcage avec champs mensuels ! lecture fichier topo correspondant a chaque snapshot climatique ! 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 06/2015 USE module3D_phy,only:nx,ny,S,slv,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp,time !use interface_input use netcdf use io_netcdf_grisli !USE printtable implicit none ! 1=decalaration variables !------------------------- integer :: nft ! NFT est le nombre de lignes a lire dans le fichier contenant le forcage climatique integer,parameter :: mois=12 integer,parameter :: ntr=1 ! nb de snapshots selon les paramètres orbitaux integer,parameter :: gtr=1 ! nb de snapshots selon l'état de la calotte integer,parameter :: ctr=1 ! nb de snapshots selon le CO2 real :: CO2_value=1120. 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,ctr,gtr) :: Ssnap ! altitude surface dans le snapshot 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 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,dimension(nx,ny,mois,ctr) :: Tm_time_fin_2 ! temperature mensuelle au temps time (corrige altitude) real,dimension(nx,ny,mois,ctr) :: Pm_time_fin_2 ! precipitation mensuelle au temps time !interpolation sur le CO2 real,dimension(nx,ny,mois) :: Tm_time_fin_3 ! temperature mensuelle au temps time (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) :: ZS !< surface topography above sea level 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 ! fichier snapshot temp file name character(len=100),dimension(ctr,gtr,ntr) :: filtr_p ! fichierprecip file CHARACTER(len=100),dimension(ctr,gtr) :: file_topo ! fichier altitude surface dans le snapshot : topo GCM sur grille GRISLI 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 :: mo,ti,tj character(len=8) :: control ! label to check clim. forc. file (filin) is usable integer :: l ! In snapshot files:the first column is the mask, read but not used real :: T_surf_ref ! variable de travail calcul temp a l'instant t a la surface S integer :: intr integer :: igtr integer :: ictr integer :: i,j character(len=100) :: file_ncdf !< fichier netcdf issue des fichiers .dat 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 DO igtr=1,gtr DO intr=1,ntr !temperature filin=TRIM(dirnameinp)//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)//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)//TRIM(file_topo(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 ! READ(20,*) Ssnap(i,j,ictr,igtr) ! if (Ssnap(i,j,ictr,igtr).eq.0.0) Ssnap(i,j,ictr,igtr)=25.0 ! END DO ! END DO ! CLOSE(20) 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_insol/filtr_t,filtr_p,file_topo ! ce bloc est a dupliquer pour chaque snapshot en changeant ! la numerotation. ntr snapshots !namelist/snap_forcage_mois/filtr_t,filtr_p,Seuil_haut,Seuil_bas,summorb,palier_ice,surf_ice,palier_CO2 ! 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_insol) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&snap_forcage_mois_insol ! 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,A)') 'file_topo = ', file_topo ! 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(f3.1,","))') 'palier_CO2 = ', palier_CO2(:) write(num_rep_42,*)'/' write(num_rep_42,428) '! fichiers temperature et precip : 12 mois et topo' write(num_rep_42,428) '! faire un bloc namelist par snapshot' 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 real COEFT,COEFP ! !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 integer :: i,j !***************** !***** ORBIT ***** !***************** Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,1) Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,1) !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(slv(i,j),S(i,j)) !Il faut mettre S0(i,j) si pas d'interpolation sinon Ssnap(i,j,ictr,igtr) !if (Zs(i,j).ge.Ssnap(i,j,ictr,igtr)) then do mo=1,mois ! correction d'altitude ! Tm_surf_mod(i,j,mo,1,1)=-lapserate(i,j,mo)*(Zs(i,j)-S0(i,j)) & ! +Tm_time_fin_1(i,j,mo,1,1) Tm_surf_mod(i,j,mo,1,1)=-lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,1,1)) & + Tm_time_fin_1(i,j,mo,1,1) ! if (Ssnap(i,j,1,1).le.25.11479) Tm_surf_mod(i,j,mo,1,1)=Tm_time_fin_1(i,j,mo,1,1) ! if (Ssnap(i,j,1,1).eq.0.) Tm_surf_mod(i,j,mo,1,1)=0.0 Pm_surf_mod(i,j,mo,1,1)=Pm_time_fin_1(i,j,mo,1,1)*exp(0.05*(Tm_surf_mod(i,j,mo,1,1) & -Tm_time_fin_1(i,j,mo,1,1))) ! if ((Ssnap(i,j,1,1).eq.0.0).and.(Zs(i,j).eq.slv(i,j))) Pm_surf_mod(i,j,mo,1,1)=0.0 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 !*************** !***** ICE ***** !*************** Tm_time_fin_2(:,:,:,:)=Tm_surf_mod(:,:,:,:,1) Pm_time_fin_2(:,:,:,:)=Pm_surf_mod(:,:,:,:,1) !*************** !***** CO2 ***** !*************** Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,1) Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,1) !*************** !**** MONTH **** !*************** 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 !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 ! calcul de Tann et Tjuly pour les sorties : Tann(:,:)=sum(Tmois,dim=3)/12. ! moy annuelle Tjuly(:,:)=Tmois(:,:,7) ! temp juillet acc(:,:)=sum(Pm_surf,dim=3,mask=Tmois < psolid) ! /12. acc(:,:)=acc(:,:)*1000./ro 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_oneway