source: branches/iLoveclim/SOURCES/climat_forcage_mod.f90 @ 359

Last change on this file since 359 was 359, checked in by aquiquet, 2 years ago

Trunk merged to branch at revision 357
File size: 17.0 KB
Line 
1module  climat_forcage_mod
2
3! Forcage climatique avec climat GCM et fichier de forcage temporel index
4! lecture de fichiers Snapshots climatiques à temps t
5! lecture d'un fichier de forcage (temp carotte de glace)
6! possibilite d'utiliser autant de snapshots que l'on veut (ntr)
7! fichiers de forcage, Snapshot et valeurs de lapse rate definis dans fichier param
8  use module3d_phy,only:nx,ny,sealevel,sealevel_2d,S,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,tafor,coefbmshelf,PI,ro,time,dirforcage,dirnameinp
9  use netcdf
10  use io_netcdf_grisli
11
12
13  implicit none
14  character(len=100) :: filin
15  integer :: nft             !   nft est le nombre de lignes a lire dans le fichier
16                           ! contenant le forcage climatique
17  real,dimension(:),allocatable :: tdate          ! time for climate forcing
18  real,dimension(:),allocatable :: tpert
19  real,dimension(:),allocatable :: spert
20  real,dimension(:),allocatable :: bpert
21 
22  integer :: ntr         ! nb of snapshot files now explicitely specified
23
24  character(len=100) :: clim_ref_file ! climat de reference
25  character(len=70), dimension(5) :: forcage_file ! fichier de forcage (climat+topo) => dimension >= ntr
26!  character(len=100) :: forcage_file2 ! fichier de forcage 2 (climat+topo)
27
28  character(len=100) :: filforc       ! nom du fichier forcage
29  real :: coefT                    !< pour modifier l'amplitude de la perturb. T
30  integer :: pertbmb               !< boolean: do we modify the bmelt under ice shelves?
31  real    :: coefbmb               !< if we modify the bmelt, we read a perturbation (based on dD, insolation, etc.) and apply this coefficient
32  real    :: r_atmvar              !< a ratio to account fast atmos variability between the 2 snapshots (in %: 0=>no fast var,1=>only fast var)
33
34
35
36!  character(len=200) :: filtr(ntr)   ! fichiers de forcage
37
38  real,dimension(:), allocatable :: ttr         !< date des tranches (snapshots)
39  real,dimension(nx,ny) :: delTatime, delTjtime, rapactime
40  real,dimension(nx,ny) :: Zs        !< surface topography above sea level
41
42  integer :: typerun ! 1 for steady state using snap 1,
43                     ! 2 for steady state using snap 2 ...
44                     ! 0 for transient
45 
46! variables pour climato mensuelle : 
47  real,dimension(nx,ny,12) :: t2m0   !< initial monthly air temperature (climato)
48  real,dimension(nx,ny,12) :: pr0    !< initial monthly precip (climato)
49  real,dimension(nx,ny) :: orog0     !< topo de la climato
50 
51! climat des snapshots :
52  real,dimension(:,:,:,:), allocatable :: t2m_ss  !< t2m monthly for every GCM snapshot
53  real,dimension(:,:,:,:), allocatable :: pr_ss   !< precip monthly for every GCM snapshot
54  real,dimension(:,:,:), allocatable :: orog_ss    !< topo for every GCM snapshot
55 
56  real,dimension(:,:,:,:), allocatable :: t2m_sstopo0 !< t2m GCM snapshot on orog0
57  real,dimension(:,:,:,:), allocatable :: pr_sstopo0  !< pr GCM snapshot on orog0
58
59
60!~   real,dimension(nx,ny,12,ntr) :: t2m_ss  !< t2m monthly for every GCM snapshot
61!~   real,dimension(nx,ny,12,ntr) :: pr_ss   !< precip monthly for every GCM snapshot
62!~   real,dimension(nx,ny,ntr) :: orog_ss    !< topo for every GCM snapshot
63 
64!~   real,dimension(nx,ny,12,ntr) :: t2m_sstopo0 !< t2m GCM snapshot on orog0
65!~   real,dimension(nx,ny,12,ntr) :: pr_sstopo0  !< pr GCM snapshot on orog0
66 
67  real :: PYY                        !< constante pour calcul de la fraction solide des precips
68  real :: psolid=2.                  !< temp limit between liquid and solid precip
69  real,dimension(nx,ny) :: FT        !< proportion of year with air temp below psolid
70  real :: lapserate                  !< gradient vertical de temp annuel et juillet
71  real :: rappact                    !< coef pour calcul du rapport accumulation
72  logical :: annual                  !< True si forcage annuel, False si forcage mensuel
73
74contains
75
76!------------------------------------------------------------------------------------------
77subroutine input_clim          !routine qui permet d'initialiser les variables climatiques
78   
79  implicit none
80  character(len=8) :: control     !< label to check clim. forc. file (filin) is usable
81  integer :: l                    !< in snapshot files:the first column is the mask, read but not used
82  integer :: err                  !< recuperation erreur
83  integer :: i,j,k,m
84  real*8, dimension(:,:), pointer   :: tab2d => null()   !< tableau 2d pointer
85  real*8, dimension(:,:,:), pointer :: tab3d => null()   !< tableau 3d pointer
86
87  deltatime(:,:)=0.
88  deltjtime(:,:)=0.
89  rapactime(:,:)=1.
90  annual=.true.  ! forcage avec fichier Tann et Tjuly
91
92! lecture du climat de reference : climato mensuelle
93  print*,'lecture du climat de reference: ',trim(DIRNAMEINP)//TRIM(clim_ref_file)
94  call Read_Ncdf_var('precip',trim(DIRNAMEINP)//TRIM(clim_ref_file),tab3d)
95  pr0(:,:,:)=tab3d(:,:,:)
96  call Read_Ncdf_var('t2m',trim(DIRNAMEINP)//TRIM(clim_ref_file),tab3d)
97  t2m0(:,:,:)=tab3d(:,:,:)
98  call Read_Ncdf_var('orog',trim(DIRNAMEINP)//TRIM(clim_ref_file),tab2d)
99  orog0(:,:)=tab2d(:,:)
100 
101!  filtr(1)=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file1)
102!  filtr(2)=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file2)
103
104! fichiers donnant l'evolution temporelle
105! -----------------------------------------
106! Pour rappel hemin40 :
107! fichier forcage : signal-cycle-hemin.dat
108! Pour anteis :
109! fichier forcage : forcmodif4cycles.dat
110  filin=trim(dirforcage)//trim(filforc) 
111
112
113! lecture des fichiers snapshots pour n'importe quel geoplace
114! -------------------------------------------------------------
115  do k=1,ntr
116     write(6,*) 'Read climate file :',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k))
117     call Read_Ncdf_var('pr',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k)),tab3d)
118     pr_ss(:,:,:,k)=tab3d(:,:,:)
119     call Read_Ncdf_var('tas',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k)),tab3d)
120     t2m_ss(:,:,:,k)=tab3d(:,:,:)
121     call Read_Ncdf_var('orog',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k)),tab2d)
122     orog_ss(:,:,k)=tab2d(:,:)
123!~      read(num_forc,*) ttr(k)
124!~      read(num_forc,*)
125!~      read(num_forc,*)
126!~      do j=1,ny
127!~         do i=1,nx
128!~            read(num_forc,*) l,rapact(i,j,k),delta(i,j,k),deltj(i,j,k)
129!~         end do
130!~      end do
131!~      close(num_forc)
132  end do
133
134
135! calcul de t2m et pr sur la topo des obs :
136  do m=1,12
137    do k=1,ntr
138      t2m_sstopo0(:,:,m,k)=t2m_ss(:,:,m,k)-lapserate*(orog0(:,:)-orog_ss(:,:,k))
139      pr_sstopo0(:,:,m,k)=pr_ss(:,:,m,k)*exp(rappact*(t2m_sstopo0(:,:,m,k)-t2m_ss(:,:,m,k)))
140    enddo
141  enddo
142
143
144! lecture des fichiers d'evolution pour tout geoplace
145! ----------------------------------------------------
146  open(num_forc,file=filin,status='old')
147!        print*,nft
148  read(num_forc,*) control,nft
149!  print*,'control :',control,nft
150! determination of file size (line nb), allocation of perturbation array
151
152  if (control.ne.'nb_lines') then
153     write(6,*) filin,'indiquer le nb de ligne en debut de fichier:'
154     write(6,*) 'le nb de lignes et le label de control nb_lines'
155     stop
156  endif
157
158  if (.not.allocated(tdate)) THEN
159     allocate(tdate(nft),stat=err)
160     if (err/=0) then
161        print *,"erreur à l'allocation du tableau tdate",err
162        stop 4
163     end if
164  end if
165
166  if (.not.allocated(tpert)) THEN
167     allocate(tpert(nft),stat=err)
168     if (err/=0) then
169        print *,"erreur à l'allocation du tableau tpert",err
170        stop 4
171     end if
172  end if
173
174  if (.not.allocated(bpert)) THEN
175     allocate(bpert(nft),stat=err)
176     if (err/=0) then
177        print *,"erreur à l'allocation du tableau tpert",err
178        stop 4
179     end if
180  end if
181
182  if (.not.allocated(spert)) THEN
183     allocate(spert(nft),stat=err)
184     if (err/=0) then
185        print *,"erreur à l'allocation du tableau spert",err
186        stop 4
187     end if
188  end if
189
190  do i=1,nft
191     if (pertbmb==1) then
192        read(num_forc,*) tdate(i),spert(i),tpert(i),bpert(i)
193     else
194        read(num_forc,*) tdate(i),spert(i),tpert(i)
195        bpert(i)=0.
196     endif
197  end do
198
199  close(num_forc)
200
201! Warning: when using climat_forcage we should in principle use a tpert which is a glacial index (ranging from about 0 to 1)
202!          here we should use the coefT to scale the signal if needed
203!          the coefbmb now is copy and paste from what we have done for Antarctica in Quiquet et al. GMD 2018
204
205  tpert(:)=tpert(:)*coefT
206  bpert(:)=max(1.+bpert(:)*coefbmb,0.01) !freezing is not allowed
207
208
209
210end subroutine input_clim
211
212
213!--------------------------------------------------------------------------------
214subroutine init_forclim
215
216  real,dimension(5) :: ttr_temp         !< date des tranches (snapshots)
217  integer :: err                  !< recuperation erreur
218  integer :: k
219 
220  namelist/clim_forcage/clim_ref_file,ntr,ttr_temp,forcage_file,typerun,lapserate,rappact,filforc,pertbmb,coeft,coefbmb,r_atmvar
221
222  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
223  read(num_param,clim_forcage)
224 
225  if (ntr.le.size(ttr_temp,dim=1)) then
226   print*,'climat_forcage : we will read ',ntr,' climate snapshots'
227   do k=ntr+1,size(ttr_temp,dim=1)
228      ttr_temp(k)= ttr_temp(ntr)
229      forcage_file(k)=' '
230   enddo
231  else
232   print*,'climat_forcage ERROR : ntr > ttr_temp size'
233   print*,'ttr_temp & forcage_file size must be >= ntr'
234   stop
235  end if
236 
237  write(num_rep_42,'(A)')'!___________________________________________________________' 
238  write(num_rep_42,clim_forcage)
239
240
241 
242  if (.not.allocated(ttr)) THEN
243   allocate(ttr(ntr),stat=err)
244   if (err/=0) then
245      print *,"Erreur à l'allocation du tableau ttr",err
246      stop 4
247   end if
248  end if
249
250  if (.not.allocated(t2m_ss)) THEN
251   allocate(t2m_ss(nx,ny,12,ntr),stat=err)
252   if (err/=0) then
253      print *,"Erreur à l'allocation du tableau t2m_ss",err
254      stop 4
255   end if
256  end if
257 
258  if (.not.allocated(pr_ss)) THEN
259   allocate(pr_ss(nx,ny,12,ntr),stat=err)
260   if (err/=0) then
261      print *,"Erreur à l'allocation du tableau pr_ss",err
262      stop 4
263   end if
264  end if
265 
266  if (.not.allocated(orog_ss)) THEN
267   allocate(orog_ss(nx,ny,ntr),stat=err)
268   if (err/=0) then
269      print *,"Erreur à l'allocation du tableau orog_ss",err
270      stop 4
271   end if
272  end if
273 
274  if (.not.allocated(t2m_sstopo0)) THEN
275   allocate(t2m_sstopo0(nx,ny,12,ntr),stat=err)
276   if (err/=0) then
277      print *,"Erreur à l'allocation du tableau t2m_sstopo0",err
278      stop 4
279   end if
280  end if
281 
282  if (.not.allocated(pr_sstopo0)) THEN
283   allocate(pr_sstopo0(nx,ny,12,ntr),stat=err)
284   if (err/=0) then
285      print *,"Erreur à l'allocation du tableau pr_sstopo0",err
286      stop 4
287   end if
288  end if
289 
290  do k=1,ntr
291   ttr(k)=ttr_temp(k)
292  enddo
293   
294
295!~   write(num_rep_42,'(A)')'!___________________________________________________________'
296!~   write(num_rep_42,'(A)') '&clim_forcage                                ! nom du bloc '
297!~   write(num_rep_42,*)
298!~   write(num_rep_42,'(A,A,A)') 'clim_ref_file = "',trim(clim_ref_file),'"'
299!~   write(num_rep_42,*) 'ttr = ', ttr
300!~ !  write(num_rep_42,'(A,A,A)') 'forcage_file= "',trim(forcage_file),'"'
301!~   write(num_rep_42,'(A,A,A)') 'forcage_file1 = "',trim(forcage_file1),'"'
302!~   write(num_rep_42,'(A,A,A)') 'forcage_file2 = "',trim(forcage_file2),'"'
303!~   write(num_rep_42,*) 'typerun = ', typerun
304!~   write(num_rep_42,*) 'lapserate = ', lapserate
305!~   write(num_rep_42,*) 'rappact = ', rappact
306!~   write(num_rep_42,'(A,A,A)') ' filforc      = "',trim(filforc),'"'
307!~   write(num_rep_42,*) 'pertbmb = ', pertbmb
308!~   write(num_rep_42,*) 'coefT   = ', coefT
309!~   write(num_rep_42,*) 'coefbmb = ', coefbmb
310!~   write(num_rep_42,*) 'r_atmvar = ', r_atmvar
311!~   write(num_rep_42,*)'/'                           
312!~   write(num_rep_42,*)
313
314  PYY=2.*PI/50.
315
316end subroutine init_forclim
317
318
319!---------------------------------------------------------------------
320!forcage climatique au cours du temps
321
322subroutine forclim
323
324!$ USE OMP_LIB
325
326  implicit none
327  real :: coeftime       !< coeftime is not used
328  real :: tafor
329  integer :: l                       !< dumm index for loops on snapsots files  l=itr,ntr-1
330  integer :: itr                     !< nb of the current snapshot file (change with time)
331  integer :: i,j,k,m
332  integer :: ift                     !< indice correspondant au pas de temps du fichier de forcage
333  real :: TEMP                       !< calcul du nbr de jour < psolid
334  real,dimension(nx,ny) :: deltaZs   ! diff temp par rapport a topo orog0
335  real,dimension(nx,ny,12) :: prmois ! precip sur topo Zs
336 
337  real (kind=kind(0.d0)) ::  timeloc
338 
339  if (typerun.eq.0) then
340     timeloc = time
341  else if ( (typerun .gt. 0) .and. (typerun .le. ubound(ttr,dim=1)) ) then
342     timeloc = ttr(typerun)
343  else
344     write(*,*) "Unknown type of simulation for climat_forcage, abort"
345     STOP
346  end if
347
348  if(timeloc.le.tdate(1)) then   ! time avant le debut du fichier forcage
349     sealevel=spert(1)
350     tafor=tpert(1)
351     coefbmshelf=bpert(1)
352     ift=1
353
354  else if (timeloc.ge.tdate(nft)) then     ! time apres la fin du fichier forcage
355     sealevel=spert(nft)
356     tafor=tpert(nft)
357     coefbmshelf=bpert(nft)
358     ift=nft
359
360  else                         ! time en general
361
362     ift = 1     
363     do i=ift,nft-1  ! interpolation entre i et i+1
364
365        if((timeloc.ge.tdate(i)).and.(timeloc.lt.tdate(i+1))) then  ! cas general
366           sealevel=spert(i)+(spert(i+1)-spert(i))*    &
367                (timeloc-tdate(i))/(tdate(i+1)-tdate(i))
368
369           tafor=tpert(i)+(tpert(i+1)-tpert(i))*      &
370                (timeloc-tdate(i))/(tdate(i+1)-tdate(i))
371
372           coefbmshelf=bpert(i)+(bpert(i+1)-bpert(i))*      &
373                (timeloc-tdate(i))/(tdate(i+1)-tdate(i))
374               
375           ift=i
376           goto 100
377        endif
378     end do
379  endif
380100 continue
381
382
383!same as for the GMD2018 paper:
384!afq -- as in Pollard and Deconto, melt is more efficient at high temperature... ->
385  if (coefbmshelf.gt.1.) then
386     coefbmshelf = 1. + ( coefbmshelf - 1. ) * 1.5 !afq.... tuning for the degla?...
387  end if
388       
389  sealevel_2d(:,:)=sealevel
390
391  if(timeloc.le.ttr(1)) then   ! time en dehors des limites du fichier forcage
392    coeftime=0.
393!~     do j=1,ny
394!~       do i=1,nx
395!~         deltatime(i,j)=delta(i,j,1)
396!~         deltjtime(i,j)=deltj(i,j,1)
397!~         rapactime(i,j)=rapact(i,j,1)
398!~       end do
399!~     end do
400    itr=1
401  else if (timeloc.ge.ttr(ntr)) then        ! mis a 0
402    coeftime=1.
403    itr=ntr
404!~     do j=1,ny
405!~       do i=1,nx
406!~         deltatime(i,j)=delta(i,j,ntr)
407!~         deltjtime(i,j)=deltj(i,j,ntr)
408!~         rapactime(i,j)=rapact(i,j,ntr)
409!~       end do
410!~     end do
411  else               !  interpolation entre l et l+1 : cas general
412     itr=1
413     do l=itr,ntr-1    ! interpolation entre i et i+1
414        if((timeloc.ge.ttr(l)).and.(timeloc.lt.ttr(l+1))) then  ! test tranche
415           coeftime= (timeloc-ttr(l))/(ttr(l+1)-ttr(l))
416!           do j=1,ny
417!              do i=1,nx
418!                 delTatime(i,j)=delta(i,j,l)+             &
419!                      (delta(i,j,l+1)-delta(i,j,l))*coeft
420!                 delTjtime(i,j)=deltj(i,j,l)+             &
421!                      (deltj(i,j,l+1)-deltj(i,j,l))*coeft
422!                 rapactime(i,j)=rapact(i,j,l)+            &
423!                      (rapact(i,j,l+1)-rapact(i,j,l))*coefp
424!              end do
425!           end do
426           itr=l
427        endif        ! fin du test sur la tranche
428     end do
429  endif   ! fin du test avant,apres,milieu
430
431  if (abs(Tafor).gt.2.) then
432     write(*,*) "A glacial index greater than 2? Really??"
433     STOP
434  endif
435  coeftime = (1-r_atmvar)*coeftime + r_atmvar*(1-Tafor)
436
437!if (geoplace(1:5).eq.'hemin') then
438
439  !$OMP PARALLEL
440!     calcul de Tjuly et et Tann
441  !$OMP WORKSHARE
442  Zs(:,:)=max(sealevel_2d(:,:),S(:,:))
443  deltaZs(:,:)= -lapserate*(Zs(:,:)-orog0(:,:)) ! difference de temperature entre Zs et orog0
444  !$OMP END WORKSHARE
445! correction d'altitude
446do m=1,12 
447  if ((itr.LT.ntr-1).and.(ntr.ge.3)) then
448! if the is more than 2 snapshots : time climate is mean between 2 snapshots and anomaly with snapshot ntr (ctrl) 
449    !$OMP WORKSHARE
450    Tmois(:,:,m)=  (t2m_sstopo0(:,:,m,itr)*(1-coeftime) + t2m_sstopo0(:,:,m,itr+1)*coeftime) - t2m_sstopo0(:,:,m,ntr) + t2m0(:,:,m) + deltaZs(:,:)
451    prmois(:,:,m)= pr0(:,:,m)*((pr_sstopo0(:,:,m,itr)*(1-coeftime) + pr_sstopo0(:,:,m,itr+1)*coeftime)/ pr_sstopo0(:,:,m,ntr)) * exp(rappact*(deltaZs(:,:)))
452    !$OMP END WORKSHARE
453
454! Temp et precip fonction de coeftime :
455  else if (itr.LT.ntr) then
456    !$OMP WORKSHARE
457    Tmois(:,:,m)=  (t2m_sstopo0(:,:,m,itr) - t2m_sstopo0(:,:,m,itr+1))*(1-coeftime) + t2m0(:,:,m) + deltaZs(:,:)
458    prmois(:,:,m)= pr0(:,:,m)*(coeftime+(1-coeftime)*(pr_sstopo0(:,:,m,itr)/pr_sstopo0(:,:,m,itr+1))) * exp(rappact*(deltaZs(:,:)))
459    !$OMP END WORKSHARE
460  else ! itr=ntr (time>tsnapmax)
461    !$OMP WORKSHARE
462    Tmois(:,:,m)=  t2m0(:,:,m) + deltaZs(:,:)
463    prmois(:,:,m)= pr0(:,:,m)*exp(rappact*(deltaZs(:,:)))
464    !$OMP END WORKSHARE
465  endif           
466enddo 
467
468
469
470  !$OMP WORKSHARE
471  Tann=sum(Tmois,dim=3)/12.
472  Tjuly=Tmois(:,:,7)
473 
474
475
476
477! calcul de l'accumulation de neige :
478  acc(:,:)=sum(prmois,dim=3,mask=Tmois < psolid)
479  acc(:,:)=acc(:,:)*1000./ro  ! conversion en glace
480  !$OMP END WORKSHARE
481  !$OMP END PARALLEL
482!debug
483!print*,'Tann, Tjuly',Tann(142,14),Tjuly(142,14) 
484!print*,'acc',acc(142,14)
485
486end subroutine forclim
487
488end module  climat_forcage_mod
Note: See TracBrowser for help on using the repository browser.