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

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

Ajout de climat-forcage-insolation_mod.f90 : Méthode JB avec interpolation sur CO2, insolation et différents états d englacement via snapshots climatique
File size: 17.0 KB
Line 
1!Interpolation NL/lineaire rayon 200/400 km + orbital forcing + CO2 decrease interpo lin ou
2!log
3! Ajout de l'interpolation verticale sous maille LMDZ -> GRISLI C. DUMAS Fev 2015
4! Lecture d'un fichier de topo associé au fichier climat
5
6module climat_forcage_insolation_mod_oneway                   
7
8! forcage avec champs mensuels
9! lecture fichier topo correspondant a chaque snapshot climatique
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 06/2015
13
14USE module3D_phy,only:nx,ny,S,slv,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp,time
15!use interface_input
16use netcdf
17use io_netcdf_grisli
18!USE printtable
19       
20implicit none
21
22! 1=decalaration variables
23!-------------------------
24
25integer :: nft          !   NFT est le nombre de lignes a lire dans le fichier contenant le forcage climatique
26
27                       
28integer,parameter :: mois=12
29integer,parameter :: ntr=1      ! nb de snapshots selon les paramètres orbitaux
30integer,parameter :: gtr=1      ! nb de snapshots selon l'état de la calotte
31integer,parameter :: ctr=1      ! nb de snapshots selon le CO2
32real :: CO2_value=1120.
33
34real,dimension(nx,ny,mois,ntr) :: Tm                    ! temperature mensuelle de chaque tranche
35real,dimension(nx,ny,mois,ntr) :: Pm                    ! precipitation mensuelle
36
37real,dimension(nx,ny,ctr,gtr) :: Ssnap                  ! altitude surface dans le snapshot
38real,dimension(nx,ny,mois,ctr,gtr,ntr) :: Tm_fin        ! tableau d'interpolation stylé
39real,dimension(nx,ny,mois,ctr,gtr,ntr) :: Pm_fin        ! tableau d'interpolation stylé
40
41!interpolation sur param orbitaux
42real,dimension(nx,ny,mois,ctr,gtr) :: Tm_time_fin_1     ! temperature mensuelle au temps time (non corrige altitude)
43real,dimension(nx,ny,mois,ctr,gtr) :: Pm_time_fin_1     ! precipitation mensuelle au temps time
44
45!interpolation sur surface glace
46real,dimension(nx,ny,mois,ctr) :: Tm_time_fin_2         ! temperature mensuelle au temps time (corrige altitude)
47real,dimension(nx,ny,mois,ctr) :: Pm_time_fin_2         ! precipitation mensuelle au temps time
48
49!interpolation sur le CO2
50real,dimension(nx,ny,mois) :: Tm_time_fin_3             ! temperature mensuelle au temps time (corrige altitude)
51real,dimension(nx,ny,mois) :: Pm_time_fin_3             ! precipitation mensuelle au temps time
52
53real,dimension(nx,ny,mois,ctr,gtr) :: Tm_surf_mod       ! correction altitude
54real,dimension(nx,ny,mois,ctr,gtr) :: Pm_surf_mod       ! correction altitude
55
56real,dimension(nx,ny,mois) :: Tm_surf                 ! surface temperature (after topo. correction)
57                                                       
58real,dimension(nx,ny,mois) :: Pm_surf                 ! surface precipitation (after topo. correction)
59
60real,dimension(nx,ny) :: ZS                             !< surface topography above sea level
61
62real,dimension(nx,ny,mois) :: lapserate                 ! lapse rate
63real :: psolid=2.                                      ! temp limit between liquid and solid precip
64
65character(len=150)  :: filin                            ! nom temporaire
66character(len=100) :: file_temporel                     ! forcage temporel
67character(len=100),dimension(ctr,gtr,ntr) :: filtr_t    ! fichier snapshot temp file name
68character(len=100),dimension(ctr,gtr,ntr) :: filtr_p    ! fichierprecip file
69CHARACTER(len=100),dimension(ctr,gtr) :: file_topo  ! fichier altitude surface dans le snapshot : topo GCM sur grille GRISLI
70
71
72real :: mincoefbmelt                                    ! butoirs pour coefbmshelf
73real :: maxcoefbmelt 
74
75
76
77contains
78
79! 2=lecture des inputs
80!--------------------
81
82subroutine input_clim                                   ! routine qui permet d'initialiser les variables climatiques
83! variables locales
84!-------------------
85
86  implicit none
87  integer :: mo,ti,tj 
88  character(len=8) :: control   ! label to check clim. forc. file (filin) is usable
89  integer :: l                  ! In snapshot files:the first column is the mask, read but not used
90  real :: T_surf_ref            ! variable de travail calcul temp a l'instant t a la surface S
91  integer :: intr
92  integer :: igtr
93  integer :: ictr
94  integer :: i,j
95  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
96  real*8, dimension(:,:,:), pointer :: data_3D => null() ! donnees lues dans le netcdf
97  real*8, dimension(:,:),pointer :: data_2D => null()    ! donnees lues dans le netcdf
98
99! lecture des fichiers snapshots pour tout geoplace
100! -------------------------------------------------
101write(6,*) 'fichiers snapshots'
102DO ictr=1,ctr
103  DO igtr=1,gtr 
104    DO intr=1,ntr
105
106      !temperature
107      filin=TRIM(dirnameinp)//TRIM(filtr_t(ictr,igtr,intr))
108      call Read_ncdf_var('t2m',trim(filin),data_3D)    ! Temperature
109      Tm_fin(:,:,:,ictr,igtr,intr)=data_3D(:,:,:) 
110!      WRITE(6,*) TRIM(filin)
111!      OPEN(20,file=TRIM(filin))
112!      DO j=1,ny
113!        DO i=1,nx
114      !    do J=ny,1,-1
115      !       do I=1,nx
116!          READ(20,*) ti, tj, (Tm_fin(i,j,mo,ictr,igtr,intr),mo=1,12)
117!        END DO
118!      END DO
119!      CLOSE(20)
120
121      !precipitation
122      filin=TRIM(dirnameinp)//TRIM(filtr_p(ictr,igtr,intr))
123      call Read_ncdf_var('precip',trim(filin),data_3D)    ! precipitation
124      Pm_fin(:,:,:,ictr,igtr,intr)=data_3D(:,:,:) 
125!      WRITE(6,*) TRIM(filin)
126!      OPEN(20,file=TRIM(filin))
127!      DO j=1,ny
128!        DO i=1,nx
129         !    do J=ny,1,-1
130         !       do I=1,nx
131!          READ(20,*) ti, tj, (Pm_fin(i,j,mo,ictr,igtr,intr),mo=1,12)
132!        END DO
133!      END DO
134!      CLOSE(20)
135    END DO
136
137       ! topo
138    filin=TRIM(dirnameinp)//TRIM(file_topo(ictr,igtr))
139    call Read_ncdf_var('TOPO',trim(filin),data_2D)    ! topo
140    Ssnap(:,:,ictr,igtr)=data_2D(:,:) 
141    where(Ssnap(:,:,ictr,igtr).eq.0.0)    ! Pour PLIOMIP niv marin=25m
142       Ssnap(:,:,ictr,igtr)=25.0
143    endwhere
144!    WRITE(6,*) TRIM(filin)
145!    OPEN(20,file=TRIM(filin))
146!    DO j=1,ny
147!      DO i=1,nx
148!        READ(20,*) Ssnap(i,j,ictr,igtr)
149!        if (Ssnap(i,j,ictr,igtr).eq.0.0) Ssnap(i,j,ictr,igtr)=25.0
150!      END DO
151!    END DO
152!    CLOSE(20)
153  ENDDO
154ENDDO
155 
156
157end subroutine input_clim
158!--------------------------------------------------------------------------------
159!subroutine input_climat_ref
160! quand on traite en absolu, pas besoin du climat de reference
161
162!end subroutine input_climat_ref
163
164
165
166SUBROUTINE init_forclim
167
168! fichiers snapshots
169  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
170  !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
171! forcage temporel
172!------------------
173  namelist/forc_temporel/file_temporel,mincoefbmelt,maxcoefbmelt
174
175! lecture par namelist
176!---------------------
177! formats pour les ecritures dans 42
178428 format(A)
179  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
180  read(num_param,snap_forcage_mois_insol)
181  write(num_rep_42,428)'!___________________________________________________________' 
182  write(num_rep_42,428) '&snap_forcage_mois_insol                       ! module climat-forcage-insolation_mod'
183  write(num_rep_42,'(A,A)')   'filtr_t = ', filtr_t
184  write(num_rep_42,'(A,A)')   'filtr_p = ', filtr_p
185  write(num_rep_42,'(A,A)')   'file_topo = ', file_topo
186!  write(num_rep_42,'(A,2(f7.1,","))')   'Seuil_haut = ', Seuil_haut(:)
187!  write(num_rep_42,'(A,2(f7.1,","))')   'Seuil_bas = ', Seuil_bas(:)
188!  write(num_rep_42,'(A,2(f5.1,","))')   'summorb = ', summorb(:)
189!  write(num_rep_42,'(A,2(f7.1,","))')   'palier_ice = ', palier_ice(:,:)
190!  write(num_rep_42,'(A,A)')   'surf_ice = ', surf_ice
191!  write(num_rep_42,'(A,2(f3.1,","))')   'palier_CO2 = ', palier_CO2(:)
192  write(num_rep_42,*)'/'                     
193  write(num_rep_42,428) '! fichiers temperature et precip : 12 mois et topo' 
194  write(num_rep_42,428) '! faire un bloc namelist par snapshot'
195  write(num_rep_42,*)
196
197!  do i=1,ntr
198! glaciaire
199!     write(filtr_t(i),'(A,i3,A)') filtr_t1(1:32),int(ttr(i)),filtr_t1(36:50)
200!     write(filtr_p(i),'(A,i3,A)') filtr_p1(1:34),int(ttr(i)),filtr_p1(38:52)
201!     write(filtr_t(i),'(A)') filtr_t
202!     write(filtr_p(i),'(A)') filtr_p
203!  enddo
204
205
206call lect_lapserate_months                               ! lit les lasperate mensuels
207                                                         ! pour une version spatialisee ecrire une autre routine
208! fichiers donnant l'evolution temporelle
209! ---------------------------- ------------
210
211rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
212read(num_param,forc_temporel)
213
214write(num_rep_42,428)'!___________________________________________________________' 
215write(num_rep_42,428) '&forc_temporel                            ! module climat_forcage_mois_mod'
216write(num_rep_42,'(A,A)')   'file_temporel =', file_temporel
217write(num_rep_42,*)     'mincoefbmelt      =', mincoefbmelt 
218write(num_rep_42,*)     'maxcoefbmelt      =', maxcoefbmelt 
219write(num_rep_42,*)'/'                     
220write(num_rep_42,428) '!fichier forcage temporel pour snapshot'
221write(num_rep_42,*)
222
223
224end subroutine init_forclim
225!---------------------------------------------------------------------
226
227!forcage climatique au cours du temps
228
229subroutine forclim
230
231implicit none
232
233real COEFT,COEFP       !
234!integer l              ! dumm index for loops on snapshots files  l=ITR,NTR-1
235!cdc integer itr            ! index of the current snapshot file (change with time)
236
237integer mo
238integer :: ictr
239integer :: igtr
240integer :: i,j
241
242!*****************
243!***** ORBIT *****
244!*****************
245Tm_time_fin_1(:,:,:,:,:)=Tm_fin(:,:,:,:,:,1)
246Pm_time_fin_1(:,:,:,:,:)=Pm_fin(:,:,:,:,:,1)
247
248
249!Correction d'altitude pour les états
250!il faut que Ssnap soit un tableau avec les différentes hauteurs de glace selon la
251!simulation (no ice, med ice, full ice et 2x, 2.5x, 3x)
252do j=1,ny
253   do i=1,nx
254      Zs(i,j)=max(slv(i,j),S(i,j))
255      !Il faut mettre S0(i,j) si pas d'interpolation sinon Ssnap(i,j,ictr,igtr)
256      !if (Zs(i,j).ge.Ssnap(i,j,ictr,igtr)) then
257      do mo=1,mois
258         ! correction d'altitude
259!         Tm_surf_mod(i,j,mo,1,1)=-lapserate(i,j,mo)*(Zs(i,j)-S0(i,j)) &
260!                                            +Tm_time_fin_1(i,j,mo,1,1)
261         Tm_surf_mod(i,j,mo,1,1)=-lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,1,1)) &
262                                  + Tm_time_fin_1(i,j,mo,1,1)
263
264      !  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)
265
266       ! if (Ssnap(i,j,1,1).eq.0.) Tm_surf_mod(i,j,mo,1,1)=0.0
267         
268         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) &
269                                            -Tm_time_fin_1(i,j,mo,1,1)))
270       
271      !  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
272     
273     end do
274            !else if (Zs(i,j).lt.Ssnap(i,j,ictr,igtr)) then 
275            !   do mo=1,mois
276            !      ! correction d'altitude avec condition pour que T ne soit
277            !      ! pas supérieure à la T de l'état sans glace corrigé
278            !      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), &
279            !                                -lapserate(i,j,mo)*(Zs(i,j)-Ssnap(i,j,ictr,gtr))+Tm_time_fin_1(i,j,mo,ictr,gtr))
280
281            !      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) &
282            !                                -Tm_time_fin_1(i,j,mo,ictr,igtr)))
283            !   end do
284            !endif
285   end do
286end do
287
288
289!***************
290!***** ICE *****
291!***************
292Tm_time_fin_2(:,:,:,:)=Tm_surf_mod(:,:,:,:,1)
293Pm_time_fin_2(:,:,:,:)=Pm_surf_mod(:,:,:,:,1)
294
295
296!***************
297!***** CO2 *****
298!***************
299Tm_time_fin_3(:,:,:)=Tm_time_fin_2(:,:,:,1)
300Pm_time_fin_3(:,:,:)=Pm_time_fin_2(:,:,:,1)
301
302
303!***************
304!**** MONTH ****
305!***************
306do i=1,nx
307   do j=1,ny
308      do mo=1,mois
309         Tm_surf(i,j,mo)=Tm_time_fin_3(i,j,mo)
310         Pm_surf(i,j,mo)=Pm_time_fin_3(i,j,mo)
311      end do
312   end do
313end do
314
315
316!coefbmshelf coefficient pour la fusion basale sous les ice shelves
317
318coefbmshelf=1.0
319coefbmshelf=max(coefbmshelf,mincoefbmelt)
320coefbmshelf=min(coefbmshelf,maxcoefbmelt)
321
322
323!*********************************
324!Correction d'altitude à vérifier
325!*********************************
326
327
328!  Correction d'altitude pour la temperature y compris sur les lacs
329! Zs est l'altitude de la surface qu'elle soit mer, glace ou lac
330
331!do j=1,ny
332!   do i=1,nx
333!      ZS(I,J)=max(slv(i,j),S(I,J))
334!      do mo=1,mois
335!         Tm_surf(i,j,mo)= - lapserate(i,j,mo) * (Zs(i,j)-Ssnap(i,j)) & ! correction d'altitude T
336!              + Tm_time(i,j,mo)
337!
338!         Pm_surf(i,j,mo)= Pm_time(i,j,mo)*exp(0.05*(Tm_surf(i,j,mo)-Tm_time(i,j,mo)))
339!
340!      end do
341!   end do
342!end do
343
344
345do mo=1,mois
346   Tmois(:,:,mo)=Tm_surf(:,:,mo)
347enddo
348
349! calcul de Tann et Tjuly pour les sorties :
350Tann(:,:)=sum(Tmois,dim=3)/12.     ! moy annuelle
351Tjuly(:,:)=Tmois(:,:,7)            ! temp juillet
352
353acc(:,:)=sum(Pm_surf,dim=3,mask=Tmois < psolid) ! /12.
354acc(:,:)=acc(:,:)*1000./ro
355
356
357END subroutine forclim
358
359
360
361!************************************************************************
362! Numerical Recipes
363! interpolation spline cubique
364
365!Récupéré grace à Christophe. Modifié 17.04.13 par JB
366SUBROUTINE splint(xa,ya,y2a,n,y)
367  INTEGER,intent(in) :: n
368  double precision,dimension(n),intent(in) :: xa,y2a,ya
369  double precision,intent(out) :: y
370! Calculates the cubic spline interpolation.
371! Given 2 arrays of dimension n, xa and ya, and y2a, the second derivative
372! of the function ya at any of the n points, it computes the interpolation for
373! the array y of dimension nmax.
374! for example, if n is 10001 (e.g. 1e6 years computed with Laskar algorithm and
375! a sampling step of 100 years), then nmax would be 1e6 because GRISLI has a
376! sampling step of 1 year.
377  INTEGER l
378  REAL a,b,h
379! We will find the right place in the table by means of bisection.
380
381! do i=1,nmax
382! print *,x(i)
383! write (*,*) 'Press Enter to Continue'
384! read (*,*)
385! enddo
386
387!l=1
388!do j=1,nmax
389!if (j==1) then
390!   l=2
391!else if (j==nmax) then
392!   l=n
393!else if (x(j).gt.xa(l)) then
394!   l=l+1
395!endif
396
397
398if (time==0) then
399   l=2
400else if (time==xa(n)) then
401   l=n
402else
403   l=1
404   do while (time.gt.xa(l))
405   l=l+1
406   enddo
407endif
408
409
410h=xa(l)-xa(l-1)
411a=(xa(l)-time)/h
412b=(time-xa(l-1))/h
413y=a*ya(l-1)+b*ya(l)+((a**3-a)*y2a(l-1)+(b**3-b)*y2a(l))*(h**2)/6.
414
415! print *,time,l,xa(l-1),xa(l),ya(l-1),ya(l),h,a,b,y
416! write (*,*) 'Press Enter to Continue'
417! read (*,*)
418
419return
420
421END SUBROUTINE splint
422
423
424
425! Calculates the 2nd derivative of the y function for any points
426! Recupere grace à Christophe. Modifié par JB 18.04.13
427SUBROUTINE spline(x,y,n,yp1,ypn,y2)
428  INTEGER,intent(in) :: n
429  double precision,dimension(n), intent(in) :: x,y
430  double precision,dimension(n), intent(out) :: y2
431  double precision,intent(in) :: yp1,ypn
432
433  INTEGER i,k
434  double precision :: p,qn,sig,un
435  double precision,dimension(n) :: u
436
437if (yp1.gt..99e30) then ! The lower boundary condition is set either to 0
438   y2(1)=0.
439   u(1)=0.
440else ! or else to have a specified first derivative.
441   y2(1)=-0.5
442   u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
443endif
444do i=2,n-1
445! This is the decomposition loop of the tridiagonal algorithm. y2 and u are used
446! for temporary storage of the decomposed factors.
447
448   sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
449   p=sig*y2(i-1)+2.
450   y2(i)=(sig-1.)/p
451   u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
452        /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
453enddo
454
455if (ypn.gt..99e30) then ! The upper boundary condition is set eith
456   qn=0.
457   un=0.
458else ! or else to have a specified first derivative.
459   qn=0.5
460   un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
461endif
462y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
463do k=n-1,1,-1 ! This is the backsubstitution loop of the tridiagonal algorithm
464   y2(k)=y2(k)*y2(k+1)+u(k)                             
465enddo
466return
467
468END SUBROUTINE spline
469
470
471
472
473
474
475
476!------------------------------------------------------------------------------------------------------------
477subroutine lect_lapserate_months  ! lapserates mensuels mais uniformes spatialement
478
479implicit none
480real,dimension(12) :: lect_lapse       ! pour la lecture
481integer :: i,j
482
483namelist/lapse_month/lect_lapse
484
485! lecture de la namelist
486
487! formats pour les ecritures dans 42
488428 format(A)
489
490rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
491read(num_param,lapse_month)
492
493write(num_rep_42,428)'!___________________________________________________________' 
494write(num_rep_42,428) '&lapse_month                                          !  module climat_forcage_mois_mod'
495write(num_rep_42,'(A,12(f0.2,","))')  'lapse_month  = ', lect_lapse
496write(num_rep_42,*)'/'                     
497write(num_rep_42,428) '! laspe rates janvier -> decembre en deg/km'
498
499
500! pour repasser en deg/m et copier dans lapserate
501
502do j=1,ny
503   do i=1,nx
504      lapserate(i,j,:)=lect_lapse(:)/1000.
505   end do
506end do
507end subroutine lect_lapserate_months
508!--------------------------------------------------------------------------------------------------------
509
510end module climat_forcage_insolation_mod_oneway
Note: See TracBrowser for help on using the repository browser.