source: trunk/SOURCES/climat-forcage-stat-mois_mod-0.1.f90 @ 23

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

initial import GRISLI trunk

File size: 14.2 KB
Line 
1!> \file climat-forcage-stat-mois_mod-0.1.f90
2!! Module pour forcage avec champs mensuels stationnaire
3!<
4
5!> \namespace climat_forcage_stat_mois_mod
6!!Module pour forcage avec champs mensuels stationnaire
7!! \author ...
8!! \date ...
9!! @note Used modules:
10!! @note   - use module3D_phy
11!! @note   - use printtable
12!<
13module climat_forcage_stat_mois_mod
14
15       USE module3D_phy
16       USE printtable
17       
18  IMPLICIT NONE
19
20! 1=decalaration variables
21!-------------------------
22
23integer :: nft 
24integer :: mo,ti,tj,la
25integer :: lake(nx,ny)                                  !< NFT est le nombre de lignes a lire dans le fichier contenant le forcage climatique
26integer,parameter :: mois=12
27integer,parameter :: NTR=1                              !< nb of snapshot files : 1 pour stationnaire - ntr is now explicitely specified in the climat module
28                       
29real,dimension(:),allocatable :: TDATE                  !< time for climate forcing
30real,dimension(:),allocatable :: alphaT
31real,dimension(:),allocatable :: alphaP
32real,dimension(:),allocatable :: SPERT
33real ttr(ntr)                                           !< la date des tranches (snapshots)  (len=ntr)
34!reak alphaTTR(ntr) ! Le alphaT de l'index glaciologiq. aux snapshots
35
36real Tm(nx,ny,mois,ntr),Pm(nx,ny,mois,ntr)
37real Tm_time(nx,ny,mois),Pm_time(nx,ny,mois)
38real Tm_surf(nx,ny,mois+1)                              !< surface temperature (after topo. correction)
39real lapserate(nx,ny,mois)
40
41character (len=100) :: filin, filin2
42character (len=100) :: filtr_t(ntr), filtr_p(ntr),filtr_t1, filtr_p1 !< snapshot file name (len=ntr)
43
44
45
46! 2=lecure des inputs
47!--------------------
48contains
49
50!> SUBROUTINE: input_clim
51!!  Routine qui permet d'initialiser les variables climatiques
52!>
53  SUBROUTINE input_clim                                         !routine qui permet d'initialiser les variables climatiques
54
55    USE module3D_phy
56
57    IMPLICIT NONE
58
59    character(len=8) :: control                         !label to check clim. forc. file (filin) is usable
60    integer :: l                                                !In snapshot files:the first column is the mask, read but not used
61
62
63    !Pm(:,:,:,:)=0.
64    !Pm_time(:,:,:)=0.
65    !Tm(:,:,:,:)=0.
66    !Tm_time(:,:,:)=0.
67    lapserate(:,:,:)=0.008
68
69
70    !====================================================HEMINOR============
71    if (geoplace.eq.'hemin40') then
72
73       !     atmospheric temperature gradient
74       TEMPGRAD=0.008
75       TEMPGRJUL=0.0065
76
77       ! 1_fichiers snapshots pour forcage
78       !----------------------------------
79       !           filtr(1)=TRIM(DIRNAMEINP)//'FORCAGES/forcage_hemin40-21k.g40'
80       !           filtr(2)=TRIM(DIRNAMEINP)//'FORCAGES/file_regions_2act.g40'
81
82       ! 2_fichiers donnant l'evolution temporelle
83       ! -----------------------------------------
84       !        filin=TRIM(DIRNAMEINP)//'signal-cycle-hemin.dat'
85
86
87       !====================================================EURASIE============
88    elseif (geoplace.eq.'euras40') then
89
90       !     atmospheric temperature gradient
91       !       TEMPGRAD=0.008
92       !       TEMPGRJUL=0.0065
93       TEMPGRAD=0.006 ! Pour Serie 2 (grad T reduit)
94       TEMPGRJUL=0.005
95
96
97       filtr_t(1)=TRIM(DIRNAMEINP)//'/ABSOLU/forc_140ka_12temper.dat'     ! 1_fichiers snapshots pour forcage
98       filtr_p(1)=TRIM(DIRNAMEINP)//'/ABSOLU/forc_140ka_12precip.dat'
99
100       filin=TRIM(DIRNAMEINP)//'signal-seal_i_inso_90ka.data'           ! 2_fichiers donnant l'evolution temporelle
101       filin=TRIM(DIRNAMEINP)//'signal-sealev_inso_26jui06.dat'
102
103       filin2=TRIM(DIRNAMEINP)//'lak_reg40km_eu40_ij'
104
105
106       !====================================================ANT20 ET 40========
107    elseif ((geoplace.eq.'anteis1').or.(geoplace.eq.'ant20km')) then
108
109
110       !     atmospheric temperature gradient
111       TEMPGRAD=0.0085
112       TEMPGRJUL=0.0085
113
114       ! 1_fichiers snapshots pour forcage
115       !----------------------------------
116       !            filtr(1)='../INPUT-DATA/FORCAGES/forcage-ant-20k.g40'
117       print*,'lecture de ',filin,filtr_p,filtr_t
118       ! 2_fichiers donnant l'evolution temporelle
119       ! ---------------------------- ------------
120       filin=TRIM(DIRNAMEINP)//'forcmodif4cycles.dat' !forcage vostok en rapport : a creer
121
122    endif !Fin du test sur geolpace
123
124    ! 3_lecure des fichiers snapsots pour tout geoplace
125    ! -------------------------------------------------
126    write(6,*) 'fichiers snapshots'
127    do k=1,ntr
128       !temperature
129       write(6,*) filtr_t(k)
130       open(20,file=filtr_t(k))
131       do j=1,ny
132          do i=1,nx
133             read(20,*) ti, tj, (Tm(i,j,mo,k),mo=1,12)
134          end do
135       end do
136       !precipitation
137       close(20)
138       close(21)
139       write(6,*) filtr_p(k)
140       open(20,file=filtr_p(k))
141       do j=1,ny
142          do i=1,nx
143             read(20,*) ti, tj, (Pm(i,j,mo,k),mo=1,12)
144          end do
145       end do
146       close(20)
147    end do
148    !       ttr(1)=-500000. !ttr est la date des tranches 3D-1.h
149    !       ttr(2)=0
150    ! 4_lecure des fichiers d'evolution pour tout geoplace
151    ! ----------------------------------------------------
152    open(20,file=filin,status='old')
153    !        print*,nft
154    read(20,*) control,nft
155    print*,'control',control,nft
156    ! Determination of file size (line nb), allocation of perturbation array
157    if (control.ne.'nb_lines') then
158       write(6,*) filin,'indiquer le nb de ligne en debut de fichier:'
159       write(6,*) 'le nb de lignes et le label de control nb_lines'
160       stop
161    endif
162
163    if (.not.allocated(tdate)) THEN
164       allocate(TDATE(nft),stat=err)
165       if (err/=0) then
166          print *,"Erreur à l'allocation du tableau TDATE",err
167          stop 4
168       end if
169    end if
170
171    if (.not.allocated(alphat)) THEN
172       allocate(alphaT(nft),stat=err)
173       if (err/=0) then
174          print *,"Erreur à l'allocation du tableau TPERT",err
175          stop 4
176       end if
177    end if
178
179    if (.not.allocated(alphap)) THEN
180       allocate(alphap(nft),stat=err)
181       if (err/=0) then
182          print *,"Erreur à l'allocation du tableau TPERT",err
183          stop 4
184       end if
185    end if
186
187    if (.not.allocated(spert)) THEN
188       allocate(spert(nft),stat=err)
189       if (err/=0) then
190          print *,"Erreur à l'allocation du tableau SPERT",err
191          stop 4
192       end if
193    end if
194
195    do I=1,NFT
196       !          read(20,*) TDATE(I),alphaT(I),alphaP(i),SPERT(I) ! forc grip
197       !          print*,i,TDATE(I),alphaT(I),alphaP(i),SPERT(I)
198       read(20,*) TDATE(I),alphaT(I),SPERT(I)           ! forc inso
199    end do
200    !RAJOUT VINCE
201    !        do k=1,ntr
202    !         do I=1,NFT
203    !           if (TDATE(I)==ttr(k)) then
204    !             alphaTTR(k)=alphaT(I)
205    !           endif
206    !         end do
207    !        end do
208
209    alphaP(:)=1.0
210    !----------------------------masque de glace pour l'Eurasie
211    do j=3,ny-2
212       do i=3,nx-2
213          mk0(i,j)=1
214          if ((i.lt.25).and.(j.gt.35)) mk0(i,j)=0
215       end do
216    end do
217    !----------------------------pour les lacs
218    !print*, 'module forcage climatique'
219    !print*, 'ajout de lakes'
220
221    !   open(50,file=filin2)
222    !       do j=1, ny
223    !          write(*,*) 'entree boucle j'
224    !          do i=1, nx
225    !       write(*,*) 'entree boucle i'
226    !       write(*,*) i
227    !            read(50,*) k, k, lake(i,j)
228    !          enddo
229    !       enddo   
230    !   close(50)
231
232    ! A CREER ou organiser a partir des codes dans ! 2_fichiers donnant l'evolution temporelle
233
234  end subroutine input_clim
235!--------------------------------------------------------------------------------
236!> SUBROUTINE: init_forclim
237!!  Routine qui permet d'initialiser les variables climatiques au cours du temps
238!>
239SUBROUTINE init_forclim
240
241!namelist/clim_forc/
242!namelist/clim_pert/rappact,retroac,rapbmshelf,mincoefbmelt,maxcoefbmelt
243
244!rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
245!read(num_param,clim_forc)
246
247! formats pour les ecritures dans 42
248428 format(A)
249
250write(num_rep_42,428)'!___________________________________________________________' 
251write(num_rep_42,428) '&clim_forc_mensuel ! pour le forcage climatique en stationnaire'
252write(num_rep_42,*)
253end subroutine init_forclim
254!---------------------------------------------------------------------
255!forcage climatique au cours du temps
256!> SUBROUTINE: forclim 
257!!  Routine permet d'effectuer le calcul climatique au cours du temps
258!>
259subroutine forclim
260
261  implicit none
262  real COEFT,COEFP,coeftime!coeftime is not used
263  INTEGER L !dumm index for loops on snapsots files  l=ITR,NTR-1
264  INTEGER ITR !nb of the current snapshot file (change with time)
265!       time en dehors des limites du fichier forcage
266INTEGER ii,jj
267              ITR=1
268!
269!run mensuels stationnaires
270! avec 1 snapshots
271do mo=1,mois
272    Tm_time(:,:,mo)=Tm(:,:,mo,1)
273    Pm_time(:,:,mo)=Pm(:,:,mo,1)
274end do
275
276!if(time.le.tdate(1)) then   ! time avant le debut du fichier forcage
277!   sealevel=spert(1)
278!   coeft=alphat(1)
279!   coefp=alphap(1)
280!   ift=1
281!
282!else if (time.ge.tdate(nft)) then     ! time apres la fin du fichier forcage
283!   sealevel=spert(nft)
284!   coeft=alphat(nft)
285!   coefp=alphap(nft)
286!   ift=nft
287!
288!else
289!
290!
291!!A modifier
292!    ift = 1      ! modifie par SC le 24/11/99
293!   do i=ift,nft-1
294!!       print*,'ds boucle1 sur snapshots',ift,TDATE(I),TDATE(I+1),time
295!          if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then
296!
297!!        interpolation entre I et I+1 : cas general
298!!           TAFOR=TPERT(I)+(TPERT(I+1)-TPERT(I))*
299!!    &            (time-TDATE(I))/(TDATE(I+1)-TDATE(I))
300!          SEALEVEL=SPERT(I)+(SPERT(I+1)-SPERT(I))*    &
301!                 (time-tdate(I))/(tdate(I+1)-tdate(I))
302!
303!          coeft=alphaT(I)+(alphaT(I+1)-alphaT(I))*      &
304!                 (time-tdate(I))/(tdate(I+1)-tdate(I))
305!
306!!          COEFP=alphaP(I)+(alphaP(I+1)-alphaP(I))*      &
307!!                 (time-TDATE(I))/(TDATE(I+1)-TDATE(I))
308!
309!                 do k=1,ntr-1 
310!                   if((time.ge.ttr(k)).and.(time.lt.ttr(k+1))) goto 105
311!
312!                 enddo
313!105  continue
314!
315!if (k==ntr) k=ntr-1
316!!Rajout Vince 26 juin 2006
317!! coefbmshelf est donné avec l'index d'insolation à 65°N
318!          coefbmshelf=coefT
319!!On retablit l'index entre
320!
321!                  if((time.ge.ttr(1)).and.(time.lt.ttr(ntr))) then
322!                      coeft  =   ( coeft - alphaTTR(k) ) / &
323!                           ( alphaTTR(k+1)-alphaTTR(k) )
324!                  else
325!                      print*,'on est pas dans les périodes forcé'
326!                      print*,'les bornes sont ttr(1) et ttr(ntr)',ttr(1),ttr(ntr)
327!                      print*,'time=',time
328!                      stop
329!                  endif
330!
331!
332!
333!! Pour les experiences eurasiennes coefp==coeft
334!          COEFP=COEFT
335!!        SEALEVEL=-50.
336!            IFT=I
337!            goto 100
338!          endif
339!
340!        end do
341!      endif
342!100   continue
343!
344!!          print*,'coeffT P',COEFT,COEFP
345!!=forcage du climat
346!!       time en dehors des limites du fichier forcage
347!        if(time.le.Ttr(1)) then   ! mis a -500 000 ans
348!          do j=1,ny
349!            do i=1,nx
350!              delTatime(i,j)=delTa(i,j,1)
351!              delTjtime(i,j)=delTj(i,j,1)
352!              rapactime(i,j)=rapact(i,j,1)
353!            end do
354!           end do
355!              ITR=1
356!        else if (time.ge.ttr(Ntr)) then        ! mis a 0
357!          do j=1,ny
358!            do i=1,nx
359!             delTatime(i,j)=delTa(i,j,ntr)
360!             delTjtime(i,j)=delTj(i,j,ntr)
361!             rapactime(i,j)=rapact(i,j,ntr)
362!            end do
363!           end do
364!           ITR=NTR
365!        else               !  interpolation entre l et l+1 : cas general
366!
367!        itr=max(itr,1)
368!       WRITE(6,*)'itr = !',itr
369!
370!
371!!parametres du fit
372!        do l=ITR,NTR-1
373!        print*,'ds boucle2 sur snapshots',l,ttr(l)
374!           if((time.ge.ttr(l)).and.(time.lt.ttr(l+1))) then  ! test tranche
375!
376!           coeftime= (time-TTR(l))/(TTR(l+1)-TTR(l))
377!           do j=1,ny
378!            do i=1,nx
379!            delTatime(i,j)=delTa(i,j,l)+             &
380!               (delTa(i,j,l+1)-delTa(i,j,l))*coefT
381!
382!            delTjtime(i,j)=delTj(i,j,l)+             &
383!               (delTj(i,j,l+1)-delTj(i,j,l))*coefT
384!
385!            rapactime(i,j)=rapact(i,j,l)+            &
386!               (rapact(i,j,l+1)-rapact(i,j,l))*coefP
387!             end do
388!           end do
389!              print*,'l= ',l,delTj(40,30,l),delTj(40,30,l+1),delTjtime(40,30)
390!           ITR=l
391!            goto 200
392!
393!            endif        ! fin du test sur la tranche
394!         end do
395!       endif   ! fin du test avant,apres,milieu
396!
397!200   continue
398print*,'l= ',l,'itr= ',itr
399
400!--------NIVEAU DES MERS ET DES LAKES
401
402if ((time.lt.-95000.).or.(time.gt.-85500.)) then ! time SANS lacs
403                levelflot=sealevel
404else
405   do j=1,ny 
406     do i=1,nx   
407!!if (ice(i,j)==1) delTjtime(i,j)=delTjtime(i,j)-3. ! correction effet albedo
408!             if ( lake(i,j)==2 ) then     !KOMI
409!               levelflot(i,j)=100.
410!             elseif ( lake(i,j)==3 ) then !SIBERIAN
411!               levelflot(i,j)=60.
412!               if ((.not.flot(i,j)).and.(flot(i+1,j)) ) print*,i,j,h(i,j),bmelt(i,j)
413!            else
414!               levelflot(i,j)=sealevel
415!            endif
416             levelflot(i,j)=-60.
417     enddo
418   enddo 
419endif
420     levelflot(:,:)=max(levelflot(:,:),bsoc(:,:))
421   nom_table='flot'
422   call printtable_l(flot,nom_table)
423   nom_table='lvflo'
424   call printtable_r(levelflot,nom_table)
425
426!   coefbmshelf(:,:)=(1.+delTatime(:,:)/7.) 
427!coefbmshelf=1. ! modif pour test de sensibilite (tof 20 avril 01)
428!coefbmshelf=coefT
429    coefbmshelf=max(coefbmshelf,0.1)
430    coefbmshelf=min(coefbmshelf,2.)
431!tafor=??
432
433print*,'time, coefbmshelf ',time, coefbmshelf
434print*,'SEA,COEFT,COEFP', SEALEVEL,COEFT,COEFP
435!rint*,'= t,sea,tafor',time,sealevel,tafor
436print*,'==========================='
437
438!     calcul de Tjuly et et Tann
439! the calcul is now done in subroutine forclim
440      do J=1,NY
441        do I=1,NX
442           ZS(I,J)=max(levelflot(i,j),S(I,J))
443            do mo=1,mois
444           Tm_surf(i,j,mo)= - lapserate(i,j,mo) * (Zs(i,j)-S0(i,j)) & ! correction d'altitude
445                      + Tm_time(i,j,mo)
446
447            end do
448        end do
449      end do
450
451!
452  Tm_surf(:,:,mois+1)=Tm_surf(:,:,1)
453
454if ((geoplace(1:5).eq.'hemin').or.(geoplace(1:5).eq.'euras')) then
455
456call accum_month() ! ds le main
457call ablation_month()
458
459elseif ((geoplace.eq.'anteis1').or.(geoplace.eq.'ant20km')) then
460! call n'est plus dans le main
461!call massb_anteis_forcage()
462!call ablation()
463else
464 print*,"partie en travaux climat forcage auter que sur heminord"
465 print*,"geoplace ",geoplace
466endif
467! call bla2
468! call bla3
469
470END subroutine forclim
471
472end module climat_forcage_stat_mois_mod
Note: See TracBrowser for help on using the repository browser.