source: trunk/SOURCES/climat-forcage_mod-0.4.f90 @ 25

Last change on this file since 25 was 4, checked in by dumas, 10 years ago

initial import GRISLI trunk

File size: 11.0 KB
Line 
1!> \file climat-forcage_mod-0.4.f90
2!!Module pour le calcule du forcage climatique
3!<
4
5!> \namespace climat_forcage_mod
6!! Module pour le calcule du forcage climatique
7!! \author ...
8!! \date ...
9!! @note Used module
10!! @note   - use module3D_phy
11!<
12module  climat_forcage_mod
13
14! declaration variables
15
16use module3d_phy
17
18implicit none
19character(len=80) :: filin
20integer :: nft             !<   nft est le nombre de lignes a lire dans le fichier
21                           !< contenant le forcage climatique
22real,dimension(:),allocatable :: tdate          !< time for climate forcing
23real,dimension(:),allocatable :: alphat
24real,dimension(:),allocatable :: alphap
25real,dimension(:),allocatable :: spert
26
27
28real :: mincoefbmelt            !< butoirs mini
29real :: maxcoefbmelt            !< butoirs maxi de coefbmshelf
30character(len=25) :: forcage_file1 !< fichier de forcage 1 (climat)
31character(len=25) :: forcage_file2 !< fichier de forcage 2
32
33
34integer,parameter :: ntr=2 !< nb of snapshot files now explicitely specified
35character(len=80) :: filtr(ntr)
36
37real ttr(ntr)!< la date des tranches (snapshots)  (len=ntr)
38real delta(nx,ny,ntr),deltj(nx,ny,ntr),rapact(nx,ny,ntr)
39real  deltatime(nx,ny),deltjtime(nx,ny),rapactime(nx,ny)
40
41contains
42
43!------------------------------------------------------------------------------------------
44
45!> SUBROUTINE: input_clim
46!! Routine qui permet d'initialiser les variables climatiques
47!>
48
49  subroutine input_clim          !routine qui permet d'initialiser les variables climatiques
50
51    implicit none
52    character(len=8) :: control     !label to check clim. forc. file (filin) is usable
53    integer l                       !in snapshot files:the first column is the mask, read but not used
54
55    deltatime(:,:)=0.
56    deltjtime(:,:)=0.
57    rapactime(:,:)=1.
58
59    if (geoplace.eq.'heminor') then
60
61       !     atmospheric temperature gradient
62       tempgrad=0.008
63       tempgrjul=0.0065
64
65
66       !pour le nord cette partie est copiee de inputfile-hemicycle.f (3juin04)
67       ! fichiers snapshots pour forcage
68       !----------------------------------
69       filtr(1)=trim(dirnameinp)//'Snapshots-GCM/forclgm-lmd5prese.g50'
70       filtr(2)=trim(dirnameinp)//'Snapshots-GCM/forclgm-lmd5futur.g50'
71
72       ! fichiers donnant l'evolution temporelle
73
74       ! this file contains tdate(i),alphat(i),alphap(i),spert(i),i=1,nft     
75       ! -------------------hemicycle ------------
76       ! lecture du fichier scalaire (sea level)
77       !        filin=trim(dirnameinp)//'heminor-forclgm.dat'
78       filin=trim(dirnameinp)//'signal-cycle-hemin.dat'
79
80    elseif (geoplace.eq.'hemin40') then
81
82       !     atmospheric temperature gradient
83       tempgrad=0.008
84       tempgrjul=0.0065
85
86       ! fichiers snapshots pour forcage
87       !----------------------------------
88       !          filtr(1)=trim(dirnameinp)//'Snapshots-GCM/forcage_hemin40-21k.g40' ! glaciaire std
89       !           filtr(1)=trim(dirnameinp)//'Snapshots-GCM/file_regions_.g40' ! glaciaire regions vincent
90       !           filtr(2)=trim(dirnameinp)//'Snapshots-GCM/file_regions_2.g40' ! pour glaciaire permanent
91       !          filtr(2)=trim(dirnameinp)//'Snapshots-GCM/forclgm-reg_oc2-00k.g40'
92       !          filtr(2)=trim(dirnameinp)//'Snapshots-GCM/no_forcage_hn40-0k.g40' ! climat act a 0ky
93       !           filtr(2)=trim(dirnameinp)//'Snapshots-GCM/file_regions_2act.g40' ! pour climat act permanent
94
95       !cdc forcage modele IPSL LSCE Cycle
96       !            filtr(1)=TRIM(DIRNAMEINP)//'Snapshots-GCM/forcage_IPSL-21k.g40'
97       !            filtr(2)=TRIM(DIRNAMEINP)//'Snapshots-GCM/no_forcage_hn40-0k.g40' !climat act a 0ka
98       filtr(1)=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file1)
99       filtr(2)=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file2)
100       print*,'fichier de forcage 1',TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file1)
101       print*,'fichier de forcage 2',TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(forcage_file2)
102
103       ! fichiers donnant l'evolution temporelle
104       ! -----------------------------------------
105       filin=trim(dirforcage)//'signal-cycle-hemin.dat'
106
107    else if ((geoplace.eq.'anteis1').or.(geoplace.eq.'ant20km')) then
108
109
110       !     atmospheric temperature gradient
111       tempgrad=0.0085
112       tempgrjul=0.0085
113
114       ! fichiers snapshots pour forcage
115       !----------------------------------
116       filtr(1)='../input-data/Snapshots-GCM/forcage-ant-20k.g40'
117       filtr(2)='../input-data/Snapshots-GCM/forcage-ant-00k.g40'
118       print*,'lecture de ',filin,filtr
119       ! fichiers donnant l'evolution temporelle
120       ! ---------------------------- ------------
121       filin=trim(dirforcage)//'forcmodif4cycles.dat' !forcage vostok en rapport : a creer
122
123    endif !fin du test sur geoplace
124
125    ! lecture des fichiers snapsots pour n'importe quel  geoplace
126    ! -------------------------------------------------------------
127    write(6,*) 'fichiers snapshots'
128    do k=1,ntr
129       write(6,*) filtr(k)
130       open(num_forc,file=filtr(k))
131       read(num_forc,*) ttr(k)
132       read(num_forc,*)
133       read(num_forc,*)
134       do j=1,ny
135          do i=1,nx
136             read(num_forc,*) l,rapact(i,j,k),delta(i,j,k),deltj(i,j,k)
137          end do
138       end do
139       close(num_forc)
140    end do
141    ttr(1)=-500000. !ttr est la date des tranches 3d-1.h
142    ttr(2)=0
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(alphat)) THEN
167       allocate(alphat(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(alphap)) THEN
175       allocate(alphap(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       read(num_forc,*) tdate(i),alphat(i),alphap(i),spert(i)
192       !          print*,tdate(i),alphat(i),alphap(i),spert(i)
193    end do
194
195    close(num_forc)
196
197
198  end subroutine input_clim
199
200
201!--------------------------------------------------------------------------------
202!> SUBROUTINE: init_forclim
203!! Initialisation des parametres du forcage climatique au cours du temps
204!>
205
206subroutine init_forclim
207
208
209namelist/clim_forcage/forcage_file1,forcage_file2,mincoefbmelt,maxcoefbmelt
210
211rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
212read(num_param,clim_forcage)
213
214
215! formats pour les ecritures dans 42
216428 format(A)
217
218write(num_rep_42,428)'!___________________________________________________________' 
219write(num_rep_42,428) '&clim_forcage                                ! nom du bloc '
220write(num_rep_42,*)
221write(num_rep_42,*) 'forcage_file1 = ', forcage_file1
222write(num_rep_42,*) 'forcage_file2 = ', forcage_file2
223write(num_rep_42,*) 'mincoefbmelt = ', mincoefbmelt
224write(num_rep_42,*) 'maxcoefbmelt = ', maxcoefbmelt
225write(num_rep_42,*)'/'                           
226write(num_rep_42,*)
227return
228end subroutine init_forclim
229
230
231!---------------------------------------------------------------------
232!> SUBROUTINE: forclim
233!! Routine pour le calcul du forcage climatique au cours du temps
234!>
235
236subroutine forclim
237
238implicit none
239real :: coeft,coefp,coeftime       !coeftime is not used
240integer :: l                       !dumm index for loops on snapsots files  l=itr,ntr-1
241integer itr                        !nb of the current snapshot file (change with time)
242
243
244if(time.le.tdate(1)) then   ! time avant le debut du fichier forcage
245   sealevel=spert(1)
246   coeft=alphat(1)
247   coefp=alphap(1)
248   ift=1
249
250else if (time.ge.tdate(nft)) then     ! time apres la fin du fichier forcage
251   sealevel=spert(nft)
252   coeft=alphat(nft)
253   coefp=alphap(nft)
254   ift=nft
255
256else                         ! time en general
257   
258   ift = 1     
259   do i=ift,nft-1  ! interpolation entre i et i+1
260
261      if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then  ! cas general
262          sealevel=spert(i)+(spert(i+1)-spert(i))*    &
263                 (time-tdate(i))/(tdate(i+1)-tdate(i))
264
265          coeft=alphat(i)+(alphat(i+1)-alphat(i))*      &
266                 (time-tdate(i))/(tdate(i+1)-tdate(i))
267
268          coefp=alphap(i)+(alphap(i+1)-alphap(i))*      &
269                 (time-tdate(i))/(tdate(i+1)-tdate(i))
270     
271!        print*,'sea,coeft,coefp', sealevel,coeft,coefp
272          ift=i
273          goto 100
274         endif
275
276      end do
277endif
278100   continue
279
280!          print*,'coefft p',coeft,coefp
281!=forcage du climat
282
283if(time.le.ttr(1)) then   ! time en dehors des limites du fichier forcage
284   do j=1,ny
285      do i=1,nx
286         deltatime(i,j)=delta(i,j,1)
287         deltjtime(i,j)=deltj(i,j,1)
288         rapactime(i,j)=rapact(i,j,1)
289      end do
290   end do
291   itr=1
292
293else if (time.ge.ttr(ntr)) then        ! mis a 0
294   do j=1,ny
295      do i=1,nx
296         deltatime(i,j)=delta(i,j,ntr)
297         deltjtime(i,j)=deltj(i,j,ntr)
298         rapactime(i,j)=rapact(i,j,ntr)
299      end do
300   end do
301   itr=ntr
302
303else               !  interpolation entre l et l+1 : cas general
304   itr=max(itr,1)
305
306   do l=itr,ntr-1    ! interpolation entre i et i+1
307     
308      if((time.ge.ttr(l)).and.(time.lt.ttr(l+1))) then  ! test tranche
309         
310         coeftime= (time-ttr(l))/(ttr(l+1)-ttr(l))
311         do j=1,ny
312            do i=1,nx
313               deltatime(i,j)=delta(i,j,l)+             &
314                    (delta(i,j,l+1)-delta(i,j,l))*coeft
315
316               deltjtime(i,j)=deltj(i,j,l)+             &
317                    (deltj(i,j,l+1)-deltj(i,j,l))*coeft
318
319               rapactime(i,j)=rapact(i,j,l)+            &
320                    (rapact(i,j,l+1)-rapact(i,j,l))*coefp
321            end do
322         end do
323         itr=l
324         goto 200
325
326      endif        ! fin du test sur la tranche
327   end do
328endif   ! fin du test avant,apres,milieu
329
330200   continue
331
332deltatime(:,:)=0.
333deltjtime(:,:)=0.
334rapactime(:,:)=1.
335
336coefbmshelf=1. ! modif pour test de sensibilite (tof 20 avril 01)
337coefbmshelf=max(coefbmshelf,mincoefbmelt)
338coefbmshelf=min(coefbmshelf,maxcoefbmelt)
339
340if (geoplace(1:5).eq.'hemin') then
341   call accum7()
342
343else if ((geoplace.eq.'anteis1').or.(geoplace.eq.'ant20km')) then
344   call massb_anteis_forcage()
345
346else
347   print*,"partie en travaux climat forcage autre que sur heminord"
348   print*,"geoplace ",geoplace
349endif
350
351
352seasea=sealevel        ! pour des sorties
353
354end subroutine forclim
355
356end module  climat_forcage_mod
Note: See TracBrowser for help on using the repository browser.