source: trunk/SOURCES/Ant16_files/climat_InitMIP_years_perturb_mod.f90 @ 135

Last change on this file since 135 was 135, checked in by aquiquet, 7 years ago

initMIP-Antarctica: setup for forward experiments
File size: 10.1 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  character(len=100) :: file_bmelt_anom    !> nom du fichier dans lequel il y a l'anomalie bmelt
72
73  integer            :: err                ! recuperation d'erreur
74  integer            :: i,j,k
75
76  !aurel for Tafor:
77  character(len=8) :: control      !label to check clim. forc. file (filin) is usable
78  character(len=80):: filin
79
80  namelist/clim_smb_T_gen/smb_file,coef_smb_unit,temp_annual_file
81  namelist/smb_anom_initMIP/file_smb_anom,massb_time
82
83428 format(A)
84  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
85  read(num_param,clim_smb_T_gen)
86
87  write(num_rep_42,428)'!________________________________________________________________' 
88  write(num_rep_42,428)'!  module climat_InitMIP_years_mod lecture climat ref          '
89  write(num_rep_42,clim_smb_T_gen)
90  write(num_rep_42,428)'! smb_file          = fichier SMB (kg/m2/an)                     '
91  write(num_rep_42,428)'! coef_smb_unit     = coef passage m glace/an  (1/910 ou 1/918)  '
92  write(num_rep_42,428)'! temp_annual_file  = Temp moy annuelle  (°C)                    '
93  write(num_rep_42,428)'!________________________________________________________________'
94
95 
96! smb : surface mass balance
97  smb_file  = trim(dirnameinp)//trim(smb_file)
98
99  call Read_Ncdf_var('smb',smb_file,tab)
100
101
102  bm(:,:)  = tab(:,:) * coef_smb_unit
103
104  acc(:,:) = 0.
105  abl(:,:) = 0.
106
107  where (bm(:,:).gt.0.)
108     acc(:,:) = bm(:,:)   ! accumulation quand positif
109  elsewhere
110     abl(:,:) = - bm(:,:) ! ablation quand negatif
111  end where
112
113
114! surface temperature  Tann
115
116  temp_annual_file = trim(dirnameinp)//trim(temp_annual_file)
117
118  call Read_Ncdf_var('Tann',temp_annual_file,tab)
119  Tann(:,:)  = tab(:,:)
120
121  ta0(:,:)   = Tann(:,:)
122  Tjuly(:,:) = Tann(:,:)
123
124
125  ! ______ Anomalies...
126 
127  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
128  read(num_param,smb_anom_initMIP)
129
130  write(num_rep_42,428)'!_______________________________________________________________________' 
131  write(num_rep_42,428)'!  module climat_InitMIP_years_mod                                      '
132  write(num_rep_42,smb_anom_initMIP)
133  write(num_rep_42,428)'! file_smb_anom     = fichier anomalie SMB de GCM                       '
134  write(num_rep_42,428)'! massb_time        = 0:fixe, 1:anomalies                               '
135  write(num_rep_42,428)'!_______________________________________________________________________' 
136
137
138  file_smb_anom  = trim(dirnameinp)//trim(file_smb_anom)
139  call Read_Ncdf_var('asmb',file_smb_anom,tab)
140  smb_anom (:,:) = Tab(:,:) !* coef_smb_unit already in m/yr
141 
142  bm_0(:,:) = bm(:,:)
143 
144  filin=trim(dirforcage)//trim(filforc)
145
146  open(num_forc,file=filin,status='old')
147 
148  read(num_forc,*) control,nft
149 
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    ! Dimensionnement des tableaux tdate, ....
159
160  if (.not.allocated(tdate)) then
161     allocate(tdate(nft),stat=err)
162     if (err/=0) then
163        print *,"erreur a l'allocation du tableau Tdate",err
164        stop 4
165     end if
166  end if
167
168  if (.not.allocated(spert)) then
169     allocate(spert(nft),stat=err)
170     if (err/=0) then
171        print *,"erreur a l'allocation du tableau Spert",err
172        stop 4
173     end if
174  end if
175 
176  if (.not.allocated(tpert)) then
177     allocate(tpert(nft),stat=err)
178     if (err/=0) then
179        print *,"erreur a l'allocation du tableau Tpert",err
180        stop 4
181     end if
182  end if
183 
184  do i=1,nft
185     read(num_forc,*) tdate(i),spert(i),tpert(i)
186  end do
187  close(num_forc)
188 
189  tpert(:)=tpert(:)*coefT
190 
191
192   
193end subroutine input_clim
194
195!--------------------------------------------------------------------------------
196!> SUBROUTINE: init_forclim
197!! Routine qui permet d'initialiser les variables climatiques au cours du temps
198!>
199subroutine init_forclim
200
201  implicit none
202  namelist/lapse_rates/T_lapse_rate
203  namelist/clim_pert_massb/coefT,filforc,pertsmb,rapsmb
204 
205  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
206  read(num_param,lapse_rates)
207
208! formats pour les ecritures dans 42
209428 format(A)
210
211  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
212  read(num_param,lapse_rates)
213
214  write(num_rep_42,428)'!________________________________________________________________' 
215  write(num_rep_42,428)'!  module climat_InitMIP_years_mod                             '
216  write(num_rep_42,lapse_rates)
217  write(num_rep_42,428)'!T_lapse_rate       = lapse rate temp annuelle                   '
218  write(num_rep_42,428)'!________________________________________________________________' 
219
220 
221  rewind(num_param) 
222  read(num_param,clim_pert_massb)
223 
224  write(num_rep_42,428)'!___________________________________________________________'
225  write(num_rep_42,428) '&clim_pert                      ! module climat_perturb_mod '
226  write(num_rep_42,*)
227  write(num_rep_42,*) 'coefT        = ', coefT
228  write(num_rep_42,'(A,A)') ' filforc      = ', filforc
229  write(num_rep_42,*) 'pertsmb        = ', pertsmb
230  write(num_rep_42,*) 'rapsmb         = ', rapsmb
231  write(num_rep_42,*)'/'
232  write(num_rep_42,*)
233
234
235! appelle la routine de lecture des smb annuels
236  call input_clim
237
238  return
239end subroutine init_forclim
240
241!--------------------------------------------------------------------------------
242!> SUBROUTINE: forclim
243!! 
244!!  Routine qui permet le calcul climatique au cours du temps
245!!  @note Au temps considere (time) attribue les scalaires
246!!  @note   - tafor : forcage en temperature
247!!  @note   - sealevel : forcage niveau des mers
248!!  @note   - coefbmelt : forcage fusion basale ice shelves
249!>
250subroutine forclim               !  au temps considere (time)
251
252  implicit none
253
254  integer i,ift
255  real             :: coefanomtime
256
257  coefanomtime = min ( real(time/50.) , 1. ) 
258  if (massb_time.eq.0) then
259     bm(:,:) = bm_0(:,:)
260  else if (massb_time.eq.1) then
261     bm(:,:) = bm_0(:,:) + coefanomtime * smb_anom(:,:)
262  end if
263 
264  if(time.lt.tdate(1)) then
265     tafor=tpert(1)
266     sealevel=spert(1)
267     ift=1
268
269  else if (time.ge.tdate(nft)) then
270     tafor=tpert(nft)
271     sealevel=spert(nft)
272     ift=nft
273           
274  else
275     do i=1,nft-1
276        if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then ! entre i et i+1 : cas general
277           tafor=tpert(i)+(tpert(i+1)-tpert(i))*       &
278                (time-tdate(i))/(tdate(i+1)-tdate(i))
279           sealevel=spert(i)+(spert(i+1)-spert(i))*    &
280                (time-tdate(i))/(tdate(i+1)-tdate(i))
281           ift=i
282           goto 100
283        endif
284     end do
285  endif
286100 continue
287
288  Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor
289  Ts(:,:)    = Tann(:,:)
290
291     ! aurel marion dufresne: we might want to decrease the SMB during glacials..?
292  if (pertsmb.eq.1) then
293     bm(:,:) = bm_0(:,:) * exp( rapsmb *(Tann(:,:)-Ta0(:,:)))
294     if (Tafor.lt.0.) then
295        where(bm(:,:).lt.0.) bm(:,:)=min(bm(:,:)-Tafor*0.05,1.) !10 degrees less give 0.5 meter more ?
296     end if
297  end if
298       
299     ! coefmshelf est un coefficient qui fait varier bmgrz et bmshelf en fonction de tafor
300  coefbmshelf=(1.+tafor/10.)       ! coefbmshelf=0 pour tafor=-10deg
301  coefbmshelf=max(coefbmshelf,0.) 
302  coefbmshelf=min(coefbmshelf,2.)
303     
304end subroutine forclim
305
306
307end module  climat_InitMIP_years_perturb_mod
Note: See TracBrowser for help on using the repository browser.