source: branches/iLoveclim/SOURCES/Ant16_files/climat_InitMIP_years_perturb_mod.f90 @ 187

Last change on this file since 187 was 187, checked in by aquiquet, 6 years ago

Grisli-iloveclim branch merged to trunk at revision 185
File size: 10.0 KB
Line 
1!> \file climat_InitMIP_mod.f90
2! forcage avec BM
3!! Module ou les variations temporelles des variables climatiques
4!! sont directement imposees
5!<
6
7!> \namespace  climat_InitMIP_mod
8!! Module ou les variations temporelles des variables climatiques
9!! sont directement imposees en suivant le protocol initMIP
10!! \author afq
11!! \date 10 oct 17
12!! @note Used modules:
13!! @note   - use module3D_phy
14!<
15module climat_InitMIP_years_perturb_mod
16
17
18use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,num_forc,dirforcage,dirnameinp,tafor,time,sealevel,coefbmshelf
19use netcdf
20use io_netcdf_grisli
21implicit none
22
23
24real :: T_lapse_rate              !< pour la temperature
25
26! anomalies: smb and bmelt
27real,dimension(nx,ny)    :: smb_anom            !> bilan de masse des anomalies indices : nx,ny
28
29
30! declaration pour le bilan de masse et le bmelt
31real,dimension(nx,ny)                :: bm_0              !> bilan de masse initial
32
33real                                 :: coef_smb_unit     ! pour corriger l'unite
34
35real,dimension(nx,ny)                :: TA0               !> Temp annuelle sur S0
36
37
38! aurel, pour climat type perturb:
39integer nft             !> nombre de lignes a lire dans le fichier forcage climatique
40real,dimension(:),allocatable :: tdate          !> time for climate forcing
41real,dimension(:),allocatable :: tpert          !> temperature for climate forcing
42real,dimension(:),allocatable :: spert          !> sea surface perturbation
43real :: coefT                    !> pour modifier l'amplitude de la perturb. T
44character(len=120) :: filforc    !> nom du fichier forcage
45integer :: pertsmb               !> boolean: do we modify the smb?
46real    :: rapsmb                !> if we modify the smb this is the equivalent of rappact
47
48
49! pour les lectures ncdf
50real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
51Character(len=10), dimension(3)      :: dimtestname       !> pour la sortie test netcdf
52integer                              :: ncidloc           !> pour la sortie test netcdf
53integer                              :: status            !> pour la sortie test netcdf
54
55integer                              :: massb_time        !< pour selectionner le type de calcul de smb
56! smb fixe ou en appliquant les anomalies
57
58contains
59
60!--------------------------------------------------------------------------------
61!> SUBROUTINE: input_clim 
62!! Routine qui permet d'initialiser les variations temporelles des variables climatiques
63!>
64!--------------------------------------------------------------------------------
65
66subroutine input_clim 
67
68  character(len=100) :: smb_file           ! fichier smb
69  character(len=100) :: temp_annual_file   ! temperature annuelles
70  character(len=100) :: file_smb_anom      !> nom du fichier dans lequel il y a l'anomalie smb
71
72  integer            :: err                ! recuperation d'erreur
73  integer            :: i
74
75  !aurel for Tafor:
76  character(len=8) :: control      !label to check clim. forc. file (filin) is usable
77  character(len=80):: filin
78
79  namelist/clim_smb_T_gen/smb_file,coef_smb_unit,temp_annual_file
80  namelist/smb_anom_initMIP/file_smb_anom,massb_time
81
82428 format(A)
83  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
84  read(num_param,clim_smb_T_gen)
85
86  write(num_rep_42,428)'!________________________________________________________________' 
87  write(num_rep_42,428)'!  module climat_InitMIP_years_mod lecture climat ref          '
88  write(num_rep_42,clim_smb_T_gen)
89  write(num_rep_42,428)'! smb_file          = fichier SMB (kg/m2/an)                     '
90  write(num_rep_42,428)'! coef_smb_unit     = coef passage m glace/an  (1/910 ou 1/918)  '
91  write(num_rep_42,428)'! temp_annual_file  = Temp moy annuelle  (°C)                    '
92  write(num_rep_42,428)'!________________________________________________________________'
93
94 
95! smb : surface mass balance
96  smb_file  = trim(dirnameinp)//trim(smb_file)
97
98  call Read_Ncdf_var('smb',smb_file,tab)
99
100
101  bm(:,:)  = tab(:,:) * coef_smb_unit
102
103  acc(:,:) = 0.
104  abl(:,:) = 0.
105
106  where (bm(:,:).gt.0.)
107     acc(:,:) = bm(:,:)   ! accumulation quand positif
108  elsewhere
109     abl(:,:) = - bm(:,:) ! ablation quand negatif
110  end where
111
112
113! surface temperature  Tann
114
115  temp_annual_file = trim(dirnameinp)//trim(temp_annual_file)
116
117  call Read_Ncdf_var('Tann',temp_annual_file,tab)
118  Tann(:,:)  = tab(:,:)
119
120  ta0(:,:)   = Tann(:,:)
121  Tjuly(:,:) = Tann(:,:)
122
123
124  ! ______ Anomalies...
125 
126  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
127  read(num_param,smb_anom_initMIP)
128
129  write(num_rep_42,428)'!_______________________________________________________________________' 
130  write(num_rep_42,428)'!  module climat_InitMIP_years_mod                                      '
131  write(num_rep_42,smb_anom_initMIP)
132  write(num_rep_42,428)'! file_smb_anom     = fichier anomalie SMB de GCM                       '
133  write(num_rep_42,428)'! massb_time        = 0:fixe, 1:anomalies                               '
134  write(num_rep_42,428)'!_______________________________________________________________________' 
135
136
137  file_smb_anom  = trim(dirnameinp)//trim(file_smb_anom)
138  call Read_Ncdf_var('asmb',file_smb_anom,tab)
139  smb_anom (:,:) = Tab(:,:) !* coef_smb_unit already in m/yr
140 
141  bm_0(:,:) = bm(:,:)
142 
143  filin=trim(dirforcage)//trim(filforc)
144
145  open(num_forc,file=filin,status='old')
146 
147  read(num_forc,*) control,nft
148 
149    ! Determination of file size (line nb), allocation of perturbation array
150 
151  if (control.ne.'nb_lines') then
152     write(6,*) filin,'indiquer le nb de ligne en debut de fichier:'
153     write(6,*) 'le nb de lignes et le label de control nb_lines'
154     stop
155  endif
156   
157    ! Dimensionnement des tableaux tdate, ....
158
159  if (.not.allocated(tdate)) then
160     allocate(tdate(nft),stat=err)
161     if (err/=0) then
162        print *,"erreur a l'allocation du tableau Tdate",err
163        stop 4
164     end if
165  end if
166
167  if (.not.allocated(spert)) then
168     allocate(spert(nft),stat=err)
169     if (err/=0) then
170        print *,"erreur a l'allocation du tableau Spert",err
171        stop 4
172     end if
173  end if
174 
175  if (.not.allocated(tpert)) then
176     allocate(tpert(nft),stat=err)
177     if (err/=0) then
178        print *,"erreur a l'allocation du tableau Tpert",err
179        stop 4
180     end if
181  end if
182 
183  do i=1,nft
184     read(num_forc,*) tdate(i),spert(i),tpert(i)
185  end do
186  close(num_forc)
187 
188  tpert(:)=tpert(:)*coefT
189 
190
191   
192end subroutine input_clim
193
194!--------------------------------------------------------------------------------
195!> SUBROUTINE: init_forclim
196!! Routine qui permet d'initialiser les variables climatiques au cours du temps
197!>
198subroutine init_forclim
199
200  implicit none
201  namelist/lapse_rates/T_lapse_rate
202  namelist/clim_pert_massb/coefT,filforc,pertsmb,rapsmb
203 
204  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
205  read(num_param,lapse_rates)
206
207! formats pour les ecritures dans 42
208428 format(A)
209
210  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
211  read(num_param,lapse_rates)
212
213  write(num_rep_42,428)'!________________________________________________________________' 
214  write(num_rep_42,428)'!  module climat_InitMIP_years_mod                             '
215  write(num_rep_42,lapse_rates)
216  write(num_rep_42,428)'!T_lapse_rate       = lapse rate temp annuelle                   '
217  write(num_rep_42,428)'!________________________________________________________________' 
218
219 
220  rewind(num_param) 
221  read(num_param,clim_pert_massb)
222 
223  write(num_rep_42,428)'!___________________________________________________________'
224  write(num_rep_42,428) '&clim_pert                      ! module climat_perturb_mod '
225  write(num_rep_42,*)
226  write(num_rep_42,*) 'coefT        = ', coefT
227  write(num_rep_42,'(A,A)') ' filforc      = ', filforc
228  write(num_rep_42,*) 'pertsmb        = ', pertsmb
229  write(num_rep_42,*) 'rapsmb         = ', rapsmb
230  write(num_rep_42,*)'/'
231  write(num_rep_42,*)
232
233
234! appelle la routine de lecture des smb annuels
235  call input_clim
236
237  return
238end subroutine init_forclim
239
240!--------------------------------------------------------------------------------
241!> SUBROUTINE: forclim
242!! 
243!!  Routine qui permet le calcul climatique au cours du temps
244!!  @note Au temps considere (time) attribue les scalaires
245!!  @note   - tafor : forcage en temperature
246!!  @note   - sealevel : forcage niveau des mers
247!!  @note   - coefbmelt : forcage fusion basale ice shelves
248!>
249subroutine forclim               !  au temps considere (time)
250
251  implicit none
252
253  integer i,ift
254  real             :: coefanomtime
255
256  coefanomtime = min ( real(time/40.) , 1. ) 
257  if (massb_time.eq.0) then
258     bm(:,:) = bm_0(:,:)
259  else if (massb_time.eq.1) then
260     bm(:,:) = bm_0(:,:) + coefanomtime * smb_anom(:,:)
261  end if
262 
263  if(time.lt.tdate(1)) then
264     tafor=tpert(1)
265     sealevel=spert(1)
266     ift=1
267
268  else if (time.ge.tdate(nft)) then
269     tafor=tpert(nft)
270     sealevel=spert(nft)
271     ift=nft
272           
273  else
274     do i=1,nft-1
275        if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then ! entre i et i+1 : cas general
276           tafor=tpert(i)+(tpert(i+1)-tpert(i))*       &
277                (time-tdate(i))/(tdate(i+1)-tdate(i))
278           sealevel=spert(i)+(spert(i+1)-spert(i))*    &
279                (time-tdate(i))/(tdate(i+1)-tdate(i))
280           ift=i
281           goto 100
282        endif
283     end do
284  endif
285100 continue
286
287  Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor
288  Ts(:,:)    = Tann(:,:)
289
290     ! aurel marion dufresne: we might want to decrease the SMB during glacials..?
291  if (pertsmb.eq.1) then
292     bm(:,:) = bm_0(:,:) * exp( rapsmb *(Tann(:,:)-Ta0(:,:)))
293     if (Tafor.lt.0.) then
294        where(bm(:,:).lt.0.) bm(:,:)=min(bm(:,:)-Tafor*0.05,1.) !10 degrees less give 0.5 meter more ?
295     end if
296  end if
297       
298     ! coefmshelf est un coefficient qui fait varier bmgrz et bmshelf en fonction de tafor
299  coefbmshelf=(1.+tafor/10.)       ! coefbmshelf=0 pour tafor=-10deg
300  coefbmshelf=max(coefbmshelf,0.) 
301  coefbmshelf=min(coefbmshelf,2.)
302     
303end subroutine forclim
304
305
306end module  climat_InitMIP_years_perturb_mod
Note: See TracBrowser for help on using the repository browser.