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

Last change on this file since 159 was 108, checked in by dumas, 7 years ago

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