source: LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90 @ 4040

Last change on this file since 4040 was 4036, checked in by crisi, 3 years ago

correction d'un probleme d'initialisation phylmdiso + suite du nettoyage dans phylmdiso
  • Property svn:executable set to *
File size: 22.3 KB
Line 
1#ifdef ISO
2! $Id: $
3
4MODULE isotopes_mod
5USE infotrac_phy, ONLY: ntraciso,niso,indnum_fn_num,ok_isotrac,use_iso, &
6&       niso_possibles
7IMPLICIT NONE
8SAVE
9
10! contient toutes les variables isotopiques et leur initialisation
11! les routines specifiquement isotopiques sont dans
12! isotopes_routines_mod pour éviter dépendance circulaire avec
13! isotopes_verif_mod.
14
15
16! indices des isotopes
17integer, save :: iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO ! indices de 1 à niso: les isos n'existant pas sont mis à 0
18!$OMP THREADPRIVATE(iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO)
19
20integer :: iso_eau_possible,iso_HDO_possible,iso_O18_possible,iso_O17_possible,iso_HTO_possible ! indices de 1 à niso_possibles: ils correspondent aux tableaux définis dans infotrac:
21! tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
22! ce sont ces indices qui doivent être utilisés avec use_iso, puisque use_iso est défini comme DIMENSION(niso_possibles)
23parameter (iso_eau_possible=1)
24parameter (iso_HDO_possible=2)
25parameter (iso_O18_possible=3)
26parameter (iso_O17_possible=4)
27parameter (iso_HTO_possible=5)
28
29integer, save :: ntracisoOR
30!$OMP THREADPRIVATE(ntracisoOR)
31
32! variables indépendantes des isotopes
33
34real, save :: pxtmelt,pxtice,pxtmin,pxtmax
35!$OMP THREADPRIVATE(pxtmelt,pxtice,pxtmin,pxtmax)
36real, save ::  tdifexp, tv0cin, thumxt1
37!$OMP THREADPRIVATE(tdifexp, tv0cin, thumxt1)
38integer, save ::  ntot
39!$OMP THREADPRIVATE(ntot)
40real, save ::  h_land_ice
41!$OMP THREADPRIVATE(h_land_ice)
42real, save ::  P_veg
43!$OMP THREADPRIVATE(P_veg)
44real, save ::  musi,lambda_sursat
45!$OMP THREADPRIVATE(lambda_sursat)
46real, save ::  Kd
47!$OMP THREADPRIVATE(Kd)
48real, save ::  rh_cste_surf_cond,T_cste_surf_cond
49!$OMP THREADPRIVATE(rh_cste_surf_cond,T_cste_surf_cond)
50
51logical, save ::   bidouille_anti_divergence
52                ! si true, rappel régulier de xteau vers q, pour éviter dérives lentes
53!$OMP THREADPRIVATE(bidouille_anti_divergence)
54logical, save ::   essai_convergence
55                ! si false, on fait rigoureusement comme dans LMDZ sans isotopes,
56                ! meme si c'est génant pour les isotopes
57!$OMP THREADPRIVATE(essai_convergence)
58integer, save ::   initialisation_iso
59                ! 0: dans fichier
60                ! 1: R=0
61                ! 2: R selon distill rayleigh
62                ! 3: R=Rsmow
63!$OMP THREADPRIVATE(initialisation_iso)
64integer, save ::  modif_SST ! 0 par defaut, 1 si on veut modifier la sst
65                ! 2 et 3: profils de SST
66!$OMP THREADPRIVATE(modif_SST)
67real, save ::  deltaTtest ! modif de la SST, uniforme.
68!$OMP THREADPRIVATE(deltaTtest)
69integer, save ::  modif_sic ! on met des trous dans glace de mer
70!$OMP THREADPRIVATE(modif_sic)
71real, save ::  deltasic ! fraction de trous minimale
72!$OMP THREADPRIVATE(deltasic)
73real, save ::  deltaTtestpoles
74!$OMP THREADPRIVATE(deltaTtestpoles)
75real, save ::  sstlatcrit
76!$OMP THREADPRIVATE(sstlatcrit)
77real, save ::  dsstlatcrit
78!$OMP THREADPRIVATE(dsstlatcrit)
79real, save ::  deltaO18_oce
80!$OMP THREADPRIVATE(deltaO18_oce)
81integer, save ::  albedo_prescrit ! 0 par defaut
82                        ! 1 si on veut garder albedo constant
83!$OMP THREADPRIVATE(albedo_prescrit)
84real, save ::  lon_min_albedo,lon_max_albedo
85!$OMP THREADPRIVATE(lon_min_albedo,lon_max_albedo)
86real, save :: lat_min_albedo,lat_max_albedo
87!$OMP THREADPRIVATE(lat_min_albedo,lat_max_albedo)
88real, save ::  deltaP_BL,tdifexp_sol
89!$OMP THREADPRIVATE(deltaP_BL,tdifexp_sol)
90integer, save ::  ruissellement_pluie,alphak_stewart
91!$OMP THREADPRIVATE(ruissellement_pluie,alphak_stewart)
92integer, save ::  calendrier_guide
93!$OMP THREADPRIVATE(calendrier_guide)
94integer, save ::  cste_surf_cond
95!$OMP THREADPRIVATE(cste_surf_cond)
96real, save ::  mixlen
97!$OMP THREADPRIVATE(mixlen)
98integer, save ::  evap_cont_cste
99!$OMP THREADPRIVATE(evap_cont_cste)
100real, save ::  deltaO18_evap_cont,d_evap_cont
101!$OMP THREADPRIVATE(deltaO18_evap_cont,d_evap_cont)
102integer, save ::  nudge_qsol,region_nudge_qsol
103!$OMP THREADPRIVATE(nudge_qsol,region_nudge_qsol)
104integer, save :: nlevmaxO17
105!$OMP THREADPRIVATE(nlevmaxO17)
106integer, save ::  no_pce
107!       real, save :: slope_limiterxy,slope_limiterz
108!$OMP THREADPRIVATE(no_pce)
109real, save ::  A_satlim
110!$OMP THREADPRIVATE(A_satlim)
111integer, save ::  ok_restrict_A_satlim,modif_ratqs
112!$OMP THREADPRIVATE(ok_restrict_A_satlim,modif_ratqs)
113real, save ::  Pcrit_ratqs,ratqsbasnew
114!$OMP THREADPRIVATE(Pcrit_ratqs,ratqsbasnew)
115real, save ::  fac_modif_evaoce
116!$OMP THREADPRIVATE(fac_modif_evaoce)
117integer, save ::  ok_bidouille_wake
118!$OMP THREADPRIVATE(ok_bidouille_wake)
119logical ::  cond_temp_env
120!$OMP THREADPRIVATE(cond_temp_env)
121
122
123! variables tableaux fn de niso
124real, ALLOCATABLE, DIMENSION(:), save :: tnat, toce, tcorr
125!$OMP THREADPRIVATE(tnat, toce, tcorr)
126real, ALLOCATABLE, DIMENSION(:), save :: tdifrel
127!$OMP THREADPRIVATE(tdifrel)
128real, ALLOCATABLE, DIMENSION(:), save :: talph1, talph2, talph3
129!$OMP THREADPRIVATE(talph1, talph2, talph3)
130real, ALLOCATABLE, DIMENSION(:), save :: talps1, talps2
131!$OMP THREADPRIVATE(talps1, talps2)
132real, ALLOCATABLE, DIMENSION(:), save :: tkcin0, tkcin1, tkcin2
133!$OMP THREADPRIVATE(tkcin0, tkcin1, tkcin2)
134real, ALLOCATABLE, DIMENSION(:), save :: alpha_liq_sol
135!$OMP THREADPRIVATE(alpha_liq_sol)
136real, ALLOCATABLE, DIMENSION(:), save :: Rdefault, Rmethox
137!$OMP THREADPRIVATE(Rdefault, Rmethox)
138character*3, ALLOCATABLE, DIMENSION(:), save :: striso
139!$OMP THREADPRIVATE(striso)
140real, save :: fac_coeff_eq17_liq, fac_coeff_eq17_ice
141!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
142
143      real ridicule ! valeur maximale pour qu'une variable de type
144                    ! rapoport de mélange puisse être considérée comme négligeable. Si
145                    ! négligeable, alors on ne verifie pas si sa compo iso esta bérrante.
146      parameter (ridicule=1e-12)     
147!      parameter (ridicule=1)
148!
149      real ridicule_rain ! valeur limite de ridicule pour les flux de pluies (rain, zrfl...)
150      parameter (ridicule_rain=1e-8) ! en kg/s <-> 1e-3mm/day
151
152      real ridicule_evap ! valeur limite de ridicule pour les evap
153      parameter (ridicule_evap=ridicule_rain*1e-2) ! en kg/s <-> 1e-3mm/day
154
155      real ridicule_qsol ! valeur limite de ridicule pour les qsol
156      parameter (ridicule_qsol=ridicule_rain) ! en kg <-> 1e-8kg
157
158      real ridicule_snow ! valeur limite de ridicule pour les snow
159      parameter (ridicule_snow=ridicule_qsol) ! en kg/s <-> 1e-8kg
160     
161        real expb_max
162        parameter (expb_max=30.0)
163
164        ! spécifique au tritium:
165       
166
167logical, save :: ok_prod_nucl_tritium ! si oui, production de tritium par essais nucleaires
168!$OMP THREADPRIVATE(ok_prod_nucl_tritium)
169        integer nessai
170        parameter (nessai=486)
171        integer, save :: day_nucl(nessai)
172!$OMP THREADPRIVATE(day_nucl)
173        integer, save :: month_nucl(nessai)
174!$OMP THREADPRIVATE(month_nucl)
175        integer, save :: year_nucl(nessai)
176!$OMP THREADPRIVATE(year_nucl)
177        real, save :: lat_nucl(nessai)
178!$OMP THREADPRIVATE(lat_nucl)
179        real, save :: lon_nucl(nessai)
180!$OMP THREADPRIVATE(lon_nucl)
181        real, save :: zmin_nucl(nessai)
182!$OMP THREADPRIVATE(zmin_nucl)
183        real, save :: zmax_nucl(nessai)
184!$OMP THREADPRIVATE(zmax_nucl)
185        real, save :: HTO_nucl(nessai)
186!$OMP THREADPRIVATE(HTO_nucl)
187
188 
189CONTAINS
190
191  SUBROUTINE iso_init()
192      use ioipsl_getin_p_mod, ONLY : getin_p
193      implicit none
194
195! -- local variables:
196
197      integer ixt
198      ! référence O18
199      real fac_enrichoce18
200      real alpha_liq_sol_O18, &
201     &     talph1_O18,talph2_O18,talph3_O18, &
202     &     talps1_O18,talps2_O18, &
203     &     tkcin0_O18,tkcin1_O18,tkcin2_O18, &
204     &     tdifrel_O18 
205     
206      ! cas de l'O17
207      real fac_kcin
208      real pente_MWL
209      integer ierr
210     
211      logical ok_nocinsat, ok_nocinsfc !sensi test
212      parameter (ok_nocinsfc=.FALSE.)  ! if T: no kinetic effect in sfc evap
213      parameter (ok_nocinsat=.FALSE.)  ! if T: no sursaturation effect for ice
214      logical Rdefault_smow
215      parameter (Rdefault_smow=.FALSE.) ! si T: Rdefault=smow; si F: nul
216      ! pour le tritium
217      integer iessai
218
219    write(*,*) 'iso_init 219: entree'
220
221! allocations mémoire
222allocate (tnat(niso))
223allocate (toce(niso))
224allocate (tcorr(niso))
225allocate (tdifrel(niso))
226allocate (talph1(niso))
227allocate (talph2(niso))
228allocate (talph3(niso))
229allocate (talps1(niso))
230allocate (talps2(niso))
231allocate (tkcin0(niso))
232allocate (tkcin1(niso))
233allocate (tkcin2(niso))
234allocate (alpha_liq_sol(niso))
235allocate (Rdefault(niso))
236allocate (Rmethox(niso))
237allocate (striso(niso))
238
239
240!--------------------------------------------------------------
241! General:
242!--------------------------------------------------------------
243
244! -- verif du nombre d'isotopes:
245      write(*,*) 'iso_init 64: niso=',niso
246
247! init de ntracisoOR: on écrasera en cas de ok_isotrac si complications avec
248! ORCHIDEE
249      ntracisoOR=ntraciso 
250             
251! -- Type of water isotopes:
252
253        iso_eau=indnum_fn_num(1)
254        iso_HDO=indnum_fn_num(2)
255        iso_O18=indnum_fn_num(3)
256        iso_O17=indnum_fn_num(4)
257        iso_HTO=indnum_fn_num(5)
258        write(*,*) 'iso_init 59: iso_eau=',iso_eau
259        write(*,*) 'iso_HDO=',iso_HDO
260        write(*,*) 'iso_O18=',iso_O18
261        write(*,*) 'iso_O17=',iso_O17
262        write(*,*) 'iso_HTO=',iso_HTO
263        write(*,*) 'iso_init 251: use_iso=',use_iso
264
265      ! initialisation
266        lambda_sursat=0.004
267        thumxt1=0.75*1.2
268        ntot=20
269        h_land_ice=20. ! à comparer aux 3000mm de snow_max
270        P_veg=1.0
271        bidouille_anti_divergence=.false.
272        essai_convergence=.false.
273        initialisation_iso=0
274        modif_sst=0
275        modif_sic=0
276        deltaTtest=0.0
277        deltasic=0.1
278        deltaTtestpoles=0.0
279        sstlatcrit=30.0
280        deltaO18_oce=0.0
281        albedo_prescrit=0
282        lon_min_albedo=-200
283        lon_max_albedo=200
284        lat_min_albedo=-100
285        lat_max_albedo=100
286        deltaP_BL=10.0
287        ruissellement_pluie=0
288        alphak_stewart=1
289        tdifexp_sol=0.67
290        calendrier_guide=0
291        cste_surf_cond=0
292mixlen=35.0       
293evap_cont_cste=0.0
294deltaO18_evap_cont=0.0
295d_evap_cont=0.0
296nudge_qsol=0
297region_nudge_qsol=1
298nlevmaxO17=50
299no_pce=0
300A_satlim=1.0
301ok_restrict_A_satlim=0
302!        slope_limiterxy=2.0
303!        slope_limiterz=2.0
304modif_ratqs=0
305Pcrit_ratqs=500.0
306ratqsbasnew=0.05
307
308fac_modif_evaoce=1.0
309ok_bidouille_wake=0
310cond_temp_env=.false.
311! si oui, la temperature de cond est celle de l'environnement,
312! pour eviter bugs quand temperature dans ascendances convs est
313! mal calculee
314ok_prod_nucl_tritium=.false.
315
316! lecture des paramètres isotopiques:
317! pour que ça marche en openMP, il faut utiliser getin_p. Car le getin ne peut
318! être appelé que par un thread à la fois, et ça pose tout un tas de problème,
319! d'où tout un tas de magouilles comme dans conf_phys_m. A terme, tout le monde
320! lira par getin_p.
321call getin_p('lambda',lambda_sursat)
322call getin_p('thumxt1',thumxt1)
323call getin_p('ntot',ntot)
324call getin_p('h_land_ice',h_land_ice)
325call getin_p('P_veg',P_veg)
326call getin_p('bidouille_anti_divergence',bidouille_anti_divergence)
327call getin_p('essai_convergence',essai_convergence)
328call getin_p('initialisation_iso',initialisation_iso)
329!if (ok_isotrac) then     
330!if (initialisation_iso.eq.0) then
331!  call getin_p('initialisation_isotrac',initialisation_isotrac)
332!endif !if (initialisation_iso.eq.0) then
333!endif !if (ok_isotrac) then     
334call getin_p('modif_sst',modif_sst)
335if (modif_sst.ge.1) then
336call getin_p('deltaTtest',deltaTtest)
337if (modif_sst.ge.2) then
338  call getin_p('deltaTtestpoles',deltaTtestpoles)
339  call getin_p('sstlatcrit',sstlatcrit)
340#ifdef ISOVERIF
341      !call iso_verif_positif(sstlatcrit,'iso_init 107')
342      if (sstlatcrit.lt.0.0) then
343        write(*,*) 'iso_init 270: sstlatcrit=',sstlatcrit
344        stop
345      endif
346#endif             
347  if (modif_sst.ge.3) then 
348      call getin_p('dsstlatcrit',dsstlatcrit)
349#ifdef ISOVERIF
350      !call iso_verif_positif(dsstlatcrit,'iso_init 110')
351      if (sstlatcrit.lt.0.0) then
352        write(*,*) 'iso_init 279: dsstlatcrit=',dsstlatcrit
353        stop
354      endif
355#endif             
356  endif !if (modif_sst.ge.3) then
357endif !if (modif_sst.ge.2) then
358endif !  if (modif_sst.ge.1) then   
359call getin_p('modif_sic',modif_sic)
360if (modif_sic.ge.1) then
361call getin_p('deltasic',deltasic)
362endif !if (modif_sic.ge.1) then
363
364call getin_p('albedo_prescrit',albedo_prescrit)
365call getin_p('lon_min_albedo',lon_min_albedo)
366call getin_p('lon_max_albedo',lon_max_albedo)
367call getin_p('lat_min_albedo',lat_min_albedo)
368call getin_p('lat_max_albedo',lat_max_albedo)
369call getin_p('deltaO18_oce',deltaO18_oce)
370call getin_p('deltaP_BL',deltaP_BL)
371call getin_p('ruissellement_pluie',ruissellement_pluie)
372call getin_p('alphak_stewart',alphak_stewart)
373call getin_p('tdifexp_sol',tdifexp_sol)
374call getin_p('calendrier_guide',calendrier_guide)
375call getin_p('cste_surf_cond',cste_surf_cond)
376call getin_p('mixlen',mixlen)
377call getin_p('evap_cont_cste',evap_cont_cste)
378call getin_p('deltaO18_evap_cont',deltaO18_evap_cont)
379call getin_p('d_evap_cont',d_evap_cont) 
380call getin_p('nudge_qsol',nudge_qsol)
381call getin_p('region_nudge_qsol',region_nudge_qsol)
382call getin_p('no_pce',no_pce)
383call getin_p('A_satlim',A_satlim)
384call getin_p('ok_restrict_A_satlim',ok_restrict_A_satlim)
385#ifdef ISOVERIF     
386!call iso_verif_positif(1.0-A_satlim,'iso_init 158')
387      if (A_satlim.gt.1.0) then
388        write(*,*) 'iso_init 315: A_satlim=',A_satlim
389        stop
390      endif
391#endif         
392!      call getin_p('slope_limiterxy',slope_limiterxy)
393!      call getin_p('slope_limiterz',slope_limiterz)
394call getin_p('modif_ratqs',modif_ratqs)
395call getin_p('Pcrit_ratqs',Pcrit_ratqs)
396call getin_p('ratqsbasnew',ratqsbasnew)
397call getin_p('fac_modif_evaoce',fac_modif_evaoce)
398call getin_p('ok_bidouille_wake',ok_bidouille_wake)
399call getin_p('cond_temp_env',cond_temp_env)
400if (use_iso(iso_HTO_possible)) then
401  ok_prod_nucl_tritium=.true.
402  call getin_p('ok_prod_nucl_tritium',ok_prod_nucl_tritium)
403endif
404
405write(*,*) 'lambda,thumxt1=',lambda_sursat,thumxt1
406write(*,*) 'bidouille_anti_divergence=',bidouille_anti_divergence
407write(*,*) 'essai_convergence=',essai_convergence
408write(*,*) 'initialisation_iso=',initialisation_iso
409write(*,*) 'modif_sst=',modif_sst
410if (modif_sst.ge.1) then
411write(*,*) 'deltaTtest=',deltaTtest
412if (modif_sst.ge.2) then 
413write(*,*) 'deltaTtestpoles,sstlatcrit=', &
414&           deltaTtestpoles,sstlatcrit
415if (modif_sst.ge.3) then   
416 write(*,*) 'dsstlatcrit=',dsstlatcrit
417endif !if (modif_sst.ge.3) then
418endif !if (modif_sst.ge.2) then
419endif !if (modif_sst.ge.1) then
420write(*,*) 'modif_sic=',modif_sic
421if (modif_sic.ge.1) then 
422write(*,*) 'deltasic=',deltasic
423endif !if (modif_sic.ge.1) then
424write(*,*) 'deltaO18_oce=',deltaO18_oce
425write(*,*) 'albedo_prescrit=',albedo_prescrit
426if (albedo_prescrit.eq.1) then
427 write(*,*) 'lon_min_albedo,lon_max_albedo=', &
428&           lon_min_albedo,lon_max_albedo
429 write(*,*) 'lat_min_albedo,lat_max_albedo=', &
430&           lat_min_albedo,lat_max_albedo
431endif !if (albedo_prescrit.eq.1) then
432write(*,*) 'deltaP_BL,ruissellement_pluie,alphak_stewart=', &
433&       deltaP_BL,ruissellement_pluie,alphak_stewart
434write(*,*) 'cste_surf_cond=',cste_surf_cond
435write(*,*) 'mixlen=',mixlen
436write(*,*) 'tdifexp_sol=',tdifexp_sol
437write(*,*) 'calendrier_guide=',calendrier_guide
438write(*,*) 'evap_cont_cste=',evap_cont_cste
439write(*,*) 'deltaO18_evap_cont,d_evap_cont=', &
440&           deltaO18_evap_cont,d_evap_cont
441write(*,*) 'nudge_qsol,region_nudge_qsol=', &
442&  nudge_qsol,region_nudge_qsol 
443write(*,*) 'nlevmaxO17=',nlevmaxO17
444write(*,*) 'no_pce=',no_pce
445write(*,*) 'A_satlim=',A_satlim
446write(*,*) 'ok_restrict_A_satlim=',ok_restrict_A_satlim
447!      write(*,*) 'slope_limiterxy=',slope_limiterxy
448!      write(*,*) 'slope_limiterz=',slope_limiterz
449write(*,*) 'modif_ratqs=',modif_ratqs
450write(*,*) 'Pcrit_ratqs=',Pcrit_ratqs
451write(*,*) 'ratqsbasnew=',ratqsbasnew
452write(*,*) 'fac_modif_evaoce=',fac_modif_evaoce
453write(*,*) 'ok_bidouille_wake=',ok_bidouille_wake
454write(*,*) 'cond_temp_env=',cond_temp_env
455write(*,*) 'ok_prod_nucl_tritium=',ok_prod_nucl_tritium
456         
457
458!--------------------------------------------------------------
459! Parameters that do not depend on the nature of water isotopes:
460!--------------------------------------------------------------
461
462! -- temperature at which ice condensate starts to form (valeur ECHAM?):
463pxtmelt=273.15
464!      pxtmelt=273.15-10.0 ! test PHASE
465
466! -- temperature at which all condensate is ice:
467pxtice=273.15-10.0
468!      pxtice=273.15-30.0 ! test PHASE
469
470! -- minimum temperature to calculate fractionation coeff
471pxtmin=273.15-120.0 ! On ne calcule qu'au dessus de -120°C
472pxtmax=273.15+60.0 ! On ne calcule qu'au dessus de +60°C
473! remarque: les coeffs ont été mesurés seulement jusq'à -40!
474
475! -- a constant for alpha_eff for equilibrium below cloud base:
476tdifexp=0.58
477tv0cin=7.0
478
479! facteurs lambda et mu dans Si=musi-lambda*T
480musi=1.0
481if (ok_nocinsat) then
482lambda_sursat = 0.0 ! no sursaturation effect
483endif           
484
485
486! diffusion dans le sol
487Kd=2.5e-9 ! m2/s   
488
489! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
490rh_cste_surf_cond=0.6
491T_cste_surf_cond=288.0
492
493!--------------------------------------------------------------
494! Parameters that depend on the nature of water isotopes:
495!--------------------------------------------------------------
496! ** constantes locales
497fac_enrichoce18=0.0005
498! on a alors tcor018=1+fac_enrichoce18
499! tcorD=1+fac_enrichoce18*8
500! tcorO17=1+fac_enrichoce18*0.528
501alpha_liq_sol_O18=1.00291 ! valeur de Lehmann & Siegenthaler, 1991,
502  ! Journal of Glaciology, vol 37, p 23
503talph1_O18=1137.
504talph2_O18=-0.4156
505talph3_O18=-2.0667E-3
506talps1_O18=11.839
507talps2_O18=-0.028244
508tkcin0_O18 = 0.006
509tkcin1_O18 = 0.000285
510tkcin2_O18 = 0.00082
511tdifrel_O18= 1./0.9723
512
513! rapport des ln(alphaeq) entre O18 et O17
514fac_coeff_eq17_liq=0.529 ! donné par Amaelle
515!      fac_coeff_eq17_ice=0.528 ! slope MWL
516fac_coeff_eq17_ice=0.529
517
518
519write(*,*) 'iso_O18,iso_HDO,iso_eau=',iso_O18,iso_HDO,iso_eau
520do 999 ixt = 1, niso
521write(*,*) 'iso_init 80: ixt=',ixt
522
523
524! -- kinetic factor for surface evaporation:
525! (cf: kcin = tkcin0                  if |V|<tv0cin
526!      kcin = tkcin1*|Vsurf| + tkcin2 if |V|>tv0cin )
527! (Rq: formula discontinuous for |V|=tv0cin... )       
528
529! -- main:
530if (ixt.eq.iso_HTO) then ! Tritium
531  tkcin0(ixt) = 0.01056
532  tkcin1(ixt) = 0.0005016
533  tkcin2(ixt) = 0.0014432
534  tnat(ixt)=0.
535  !toce(ixt)=2.2222E-8 ! corrigé par Alex Cauquoin
536  !toce(ixt)=1.0E-18 ! rapport 3H/1H ocean
537  toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
538  tcorr(ixt)=1.
539  tdifrel(ixt)=1./0.968
540  talph1(ixt)=46480.
541  talph2(ixt)=-103.87
542  talph3(ixt)=0.
543  talps1(ixt)=46480.
544  talps2(ixt)=-103.87
545  alpha_liq_sol(ixt)=1.
546  Rdefault(ixt)=0.0
547  Rmethox(ixt)=0.0
548  striso(ixt)='HTO'
549endif     
550if (ixt.eq.iso_O17) then ! Deuterium
551  pente_MWL=0.528
552!          tdifrel(ixt)=1./0.985452 ! donné par Amaelle
553  tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG
554!          fac_kcin=0.5145 ! donné par Amaelle
555  fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
556  tkcin0(ixt) = tkcin0_O18*fac_kcin
557  tkcin1(ixt) = tkcin1_O18*fac_kcin
558  tkcin2(ixt) = tkcin2_O18*fac_kcin
559  tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
560  toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
561  tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle           
562  talph1(ixt)=talph1_O18
563  talph2(ixt)=talph2_O18
564  talph3(ixt)=talph3_O18
565  talps1(ixt)=talps1_O18
566  talps2(ixt)=talps2_O18     
567  alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq
568  if (Rdefault_smow) then   
569        Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
570  else
571        Rdefault(ixt)=0.0
572  endif
573  Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
574  striso(ixt)='O17'
575endif
576
577if (ixt.eq.iso_O18) then ! Oxygene18
578  tkcin0(ixt) = tkcin0_O18
579  tkcin1(ixt) = tkcin1_O18
580  tkcin2(ixt) = tkcin2_O18
581  tnat(ixt)=2005.2E-6
582  toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
583  tcorr(ixt)=1.0+fac_enrichoce18
584  tdifrel(ixt)=tdifrel_O18
585  talph1(ixt)=talph1_O18
586  talph2(ixt)=talph2_O18
587  talph3(ixt)=talph3_O18
588  talps1(ixt)=talps1_O18
589  talps2(ixt)=talps2_O18
590  alpha_liq_sol(ixt)=alpha_liq_sol_O18   
591  if (Rdefault_smow) then   
592        Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
593  else
594        Rdefault(ixt)=0.0
595  endif
596  Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006   
597!       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
598  striso(ixt)='O18'
599  write(*,*) 'isotopes_mod 519: ixt,striso(ixt)=',ixt,striso(ixt)
600endif
601
602if (ixt.eq.iso_HDO) then ! Deuterium
603  pente_MWL=8.0
604!          fac_kcin=0.88
605  tdifrel(ixt)=1./0.9755
606  fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
607  tkcin0(ixt) = tkcin0_O18*fac_kcin
608  tkcin1(ixt) = tkcin1_O18*fac_kcin
609  tkcin2(ixt) = tkcin2_O18*fac_kcin
610  tnat(ixt)=155.76E-6
611  toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
612  tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL         
613  talph1(ixt)=24844.
614  talph2(ixt)=-76.248
615  talph3(ixt)=52.612E-3
616  talps1(ixt)=16288.
617  talps2(ixt)=-0.0934
618  !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
619  alpha_liq_sol(ixt)=1.0212
620  ! valeur de Lehmann & Siegenthaler, 1991, Journal of
621  ! Glaciology, vol 37, p 23
622  if (Rdefault_smow) then   
623    Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
624  else
625    Rdefault(ixt)=0.0
626  endif
627  Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
628  striso(ixt)='HDO'
629  write(*,*) 'isotopes_mod 548: ixt,striso(ixt)=',ixt,striso(ixt)
630endif
631
632!       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
633if (ixt.eq.iso_eau) then ! Oxygene16
634  tkcin0(ixt) = 0.0
635  tkcin1(ixt) = 0.0
636  tkcin2(ixt) = 0.0
637  tnat(ixt)=1.
638  toce(ixt)=tnat(ixt)
639  tcorr(ixt)=1.0
640  tdifrel(ixt)=1.
641  talph1(ixt)=0.
642  talph2(ixt)=0.
643  talph3(ixt)=0.
644  talps1(ixt)=0.
645  talph3(ixt)=0.
646  alpha_liq_sol(ixt)=1.
647  if (Rdefault_smow) then
648        Rdefault(ixt)=tnat(ixt)*1.0
649  else
650        Rdefault(ixt)=1.0
651  endif
652  Rmethox(ixt)=1.0
653  striso(ixt)='eau'
654endif
655
656999   continue
657
658! test de sensibilité:
659if (ok_nocinsfc) then ! no kinetic effect in sfc evaporation
660 do ixt=1,niso
661  tkcin0(ixt) = 0.0
662  tkcin1(ixt) = 0.0
663  tkcin2(ixt) = 0.0
664 enddo
665endif
666
667! nom des isotopes
668
669! verif
670write(*,*) 'iso_init 285: verif initialisation:'
671
672do ixt=1,niso
673  write(*,*) '* striso(',ixt,')=<'//striso(ixt)//'>'
674  write(*,*) 'tnat(',ixt,')=',tnat(ixt)
675!          write(*,*) 'alpha_liq_sol(',ixt,')=',alpha_liq_sol(ixt)
676!          write(*,*) 'tkcin0(',ixt,')=',tkcin0(ixt)
677!          write(*,*) 'tdifrel(',ixt,')=',tdifrel(ixt)
678enddo
679write(*,*) 'iso_init 69: lambda=',lambda_sursat
680write(*,*) 'iso_init 69: thumxt1=',thumxt1
681write(*,*) 'iso_init 69: h_land_ice=',h_land_ice
682write(*,*) 'iso_init 69: P_veg=',P_veg   
683
684    return
685END SUBROUTINE iso_init
686
687
688END MODULE isotopes_mod
689#endif
690
691
Note: See TracBrowser for help on using the repository browser.