source: trunk/SOURCES/climat_forcage_mod.f90 @ 250

Last change on this file since 250 was 250, checked in by aquiquet, 5 years ago

Fast atmospheric variability taken into account in climat_forcage using a glacial index

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