source: LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90_bak7fev2018 @ 3927

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR
  • Property svn:executable set to *
File size: 22.4 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
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! allocations mémoire
220
221allocate (tnat(niso))
222allocate (toce(niso))
223allocate (tcorr(niso))
224allocate (tdifrel(niso))
225allocate (talph1(niso))
226allocate (talph2(niso))
227allocate (talph3(niso))
228allocate (talps1(niso))
229allocate (talps2(niso))
230allocate (tkcin0(niso))
231allocate (tkcin1(niso))
232allocate (tkcin2(niso))
233allocate (alpha_liq_sol(niso))
234allocate (Rdefault(niso))
235allocate (Rmethox(niso))
236allocate (striso(niso))
237
238
239!--------------------------------------------------------------
240! General:
241!--------------------------------------------------------------
242
243! -- verif du nombre d'isotopes:
244      write(*,*) 'iso_init 64: niso=',niso
245
246! init de ntracisoOR: on écrasera en cas de ok_isotrac si complications avec
247! ORCHIDEE
248      ntracisoOR=ntraciso 
249             
250! -- Type of water isotopes:
251
252        iso_eau=indnum_fn_num(1)
253        iso_HDO=indnum_fn_num(2)
254        iso_O18=indnum_fn_num(3)
255        iso_O17=indnum_fn_num(4)
256        iso_HTO=indnum_fn_num(5)
257        write(*,*) 'iso_init 59: iso_eau=',iso_eau
258        write(*,*) 'iso_HDO=',iso_HDO
259        write(*,*) 'iso_O18=',iso_O18
260        write(*,*) 'iso_O17=',iso_O17
261        write(*,*) 'iso_HTO=',iso_HTO
262        write(*,*) 'iso_init 251: use_iso=',use_iso
263
264      ! initialisation
265        lambda_sursat=0.004
266        thumxt1=0.75*1.2
267        ntot=20
268        h_land_ice=20. ! à comparer aux 3000mm de snow_max
269        P_veg=1.0
270        bidouille_anti_divergence=.false.
271        essai_convergence=.false.
272        initialisation_iso=0
273        modif_sst=0
274        modif_sic=0
275        deltaTtest=0.0
276        deltasic=0.1
277        deltaTtestpoles=0.0
278        sstlatcrit=30.0
279        deltaO18_oce=0.0
280        albedo_prescrit=0
281        lon_min_albedo=-200
282        lon_max_albedo=200
283        lat_min_albedo=-100
284        lat_max_albedo=100
285        deltaP_BL=10.0
286        ruissellement_pluie=0
287        alphak_stewart=1
288        tdifexp_sol=0.67
289        calendrier_guide=0
290        cste_surf_cond=0
291mixlen=35.0       
292evap_cont_cste=0.0
293deltaO18_evap_cont=0.0
294d_evap_cont=0.0
295nudge_qsol=0
296region_nudge_qsol=1
297nlevmaxO17=50
298no_pce=0
299A_satlim=1.0
300ok_restrict_A_satlim=0
301!        slope_limiterxy=2.0
302!        slope_limiterz=2.0
303modif_ratqs=0
304Pcrit_ratqs=500.0
305ratqsbasnew=0.05
306
307fac_modif_evaoce=1.0
308ok_bidouille_wake=0
309cond_temp_env=.false.
310! si oui, la temperature de cond est celle de l'environnement,
311! pour eviter bugs quand temperature dans ascendances convs est
312! mal calculee
313ok_prod_nucl_tritium=.false.
314
315! lecture des paramètres isotopiques:
316call getin('lambda',lambda_sursat)
317call getin('thumxt1',thumxt1)
318call getin('ntot',ntot)
319call getin('h_land_ice',h_land_ice)
320call getin('P_veg',P_veg)
321call getin('bidouille_anti_divergence',bidouille_anti_divergence)
322call getin('essai_convergence',essai_convergence)
323call getin('initialisation_iso',initialisation_iso)
324!if (ok_isotrac) then     
325!if (initialisation_iso.eq.0) then
326!  call getin('initialisation_isotrac',initialisation_isotrac)
327!endif !if (initialisation_iso.eq.0) then
328!endif !if (ok_isotrac) then     
329call getin('modif_sst',modif_sst)
330if (modif_sst.ge.1) then
331call getin('deltaTtest',deltaTtest)
332if (modif_sst.ge.2) then
333  call getin('deltaTtestpoles',deltaTtestpoles)
334  call getin('sstlatcrit',sstlatcrit)
335#ifdef ISOVERIF
336      !call iso_verif_positif(sstlatcrit,'iso_init 107')
337      if (sstlatcrit.lt.0.0) then
338        write(*,*) 'iso_init 270: sstlatcrit=',sstlatcrit
339        stop
340      endif
341#endif             
342  if (modif_sst.ge.3) then 
343      call getin('dsstlatcrit',dsstlatcrit)
344#ifdef ISOVERIF
345      !call iso_verif_positif(dsstlatcrit,'iso_init 110')
346      if (sstlatcrit.lt.0.0) then
347        write(*,*) 'iso_init 279: dsstlatcrit=',dsstlatcrit
348        stop
349      endif
350#endif             
351  endif !if (modif_sst.ge.3) then
352endif !if (modif_sst.ge.2) then
353endif !  if (modif_sst.ge.1) then   
354call getin('modif_sic',modif_sic)
355if (modif_sic.ge.1) then
356call getin('deltasic',deltasic)
357endif !if (modif_sic.ge.1) then
358
359call getin('albedo_prescrit',albedo_prescrit)
360call getin('lon_min_albedo',lon_min_albedo)
361call getin('lon_max_albedo',lon_max_albedo)
362call getin('lat_min_albedo',lat_min_albedo)
363call getin('lat_max_albedo',lat_max_albedo)
364call getin('deltaO18_oce',deltaO18_oce)
365call getin('deltaP_BL',deltaP_BL)
366call getin('ruissellement_pluie',ruissellement_pluie)
367call getin('alphak_stewart',alphak_stewart)
368call getin('tdifexp_sol',tdifexp_sol)
369call getin('calendrier_guide',calendrier_guide)
370call getin('cste_surf_cond',cste_surf_cond)
371call getin('mixlen',mixlen)
372call getin('evap_cont_cste',evap_cont_cste)
373call getin('deltaO18_evap_cont',deltaO18_evap_cont)
374call getin('d_evap_cont',d_evap_cont) 
375call getin('nudge_qsol',nudge_qsol)
376call getin('region_nudge_qsol',region_nudge_qsol)
377call getin('no_pce',no_pce)
378call getin('A_satlim',A_satlim)
379call getin('ok_restrict_A_satlim',ok_restrict_A_satlim)
380#ifdef ISOVERIF     
381!call iso_verif_positif(1.0-A_satlim,'iso_init 158')
382      if (A_satlim.gt.1.0) then
383        write(*,*) 'iso_init 315: A_satlim=',A_satlim
384        stop
385      endif
386#endif         
387!      call getin('slope_limiterxy',slope_limiterxy)
388!      call getin('slope_limiterz',slope_limiterz)
389call getin('modif_ratqs',modif_ratqs)
390call getin('Pcrit_ratqs',Pcrit_ratqs)
391call getin('ratqsbasnew',ratqsbasnew)
392call getin('fac_modif_evaoce',fac_modif_evaoce)
393call getin('ok_bidouille_wake',ok_bidouille_wake)
394call getin('cond_temp_env',cond_temp_env)
395if (use_iso(iso_HTO_possible)) then
396  ok_prod_nucl_tritium=.true.
397  call getin('ok_prod_nucl_tritium',ok_prod_nucl_tritium)
398endif
399
400write(*,*) 'lambda,thumxt1=',lambda_sursat,thumxt1
401write(*,*) 'bidouille_anti_divergence=',bidouille_anti_divergence
402write(*,*) 'essai_convergence=',essai_convergence
403write(*,*) 'initialisation_iso=',initialisation_iso
404write(*,*) 'modif_sst=',modif_sst
405if (modif_sst.ge.1) then
406write(*,*) 'deltaTtest=',deltaTtest
407if (modif_sst.ge.2) then 
408write(*,*) 'deltaTtestpoles,sstlatcrit=', &
409&           deltaTtestpoles,sstlatcrit
410if (modif_sst.ge.3) then   
411 write(*,*) 'dsstlatcrit=',dsstlatcrit
412endif !if (modif_sst.ge.3) then
413endif !if (modif_sst.ge.2) then
414endif !if (modif_sst.ge.1) then
415write(*,*) 'modif_sic=',modif_sic
416if (modif_sic.ge.1) then 
417write(*,*) 'deltasic=',deltasic
418endif !if (modif_sic.ge.1) then
419write(*,*) 'deltaO18_oce=',deltaO18_oce
420write(*,*) 'albedo_prescrit=',albedo_prescrit
421if (albedo_prescrit.eq.1) then
422 write(*,*) 'lon_min_albedo,lon_max_albedo=', &
423&           lon_min_albedo,lon_max_albedo
424 write(*,*) 'lat_min_albedo,lat_max_albedo=', &
425&           lat_min_albedo,lat_max_albedo
426endif !if (albedo_prescrit.eq.1) then
427write(*,*) 'deltaP_BL,ruissellement_pluie,alphak_stewart=', &
428&       deltaP_BL,ruissellement_pluie,alphak_stewart
429write(*,*) 'cste_surf_cond=',cste_surf_cond
430write(*,*) 'mixlen=',mixlen
431write(*,*) 'tdifexp_sol=',tdifexp_sol
432write(*,*) 'calendrier_guide=',calendrier_guide
433write(*,*) 'evap_cont_cste=',evap_cont_cste
434write(*,*) 'deltaO18_evap_cont,d_evap_cont=', &
435&           deltaO18_evap_cont,d_evap_cont
436write(*,*) 'nudge_qsol,region_nudge_qsol=', &
437&  nudge_qsol,region_nudge_qsol 
438write(*,*) 'nlevmaxO17=',nlevmaxO17
439write(*,*) 'no_pce=',no_pce
440write(*,*) 'A_satlim=',A_satlim
441write(*,*) 'ok_restrict_A_satlim=',ok_restrict_A_satlim
442!      write(*,*) 'slope_limiterxy=',slope_limiterxy
443!      write(*,*) 'slope_limiterz=',slope_limiterz
444write(*,*) 'modif_ratqs=',modif_ratqs
445write(*,*) 'Pcrit_ratqs=',Pcrit_ratqs
446write(*,*) 'ratqsbasnew=',ratqsbasnew
447write(*,*) 'fac_modif_evaoce=',fac_modif_evaoce
448write(*,*) 'ok_bidouille_wake=',ok_bidouille_wake
449write(*,*) 'cond_temp_env=',cond_temp_env
450write(*,*) 'ok_prod_nucl_tritium=',ok_prod_nucl_tritium
451         
452
453!--------------------------------------------------------------
454! Parameters that do not depend on the nature of water isotopes:
455!--------------------------------------------------------------
456
457! -- temperature at which ice condensate starts to form (valeur ECHAM?):
458pxtmelt=273.15
459!      pxtmelt=273.15-10.0 ! test PHASE
460
461! -- temperature at which all condensate is ice:
462pxtice=273.15-10.0
463!      pxtice=273.15-30.0 ! test PHASE
464
465! -- minimum temperature to calculate fractionation coeff
466pxtmin=273.15-120.0 ! On ne calcule qu'au dessus de -120°C
467pxtmax=273.15+60.0 ! On ne calcule qu'au dessus de +60°C
468! remarque: les coeffs ont été mesurés seulement jusq'à -40!
469
470! -- a constant for alpha_eff for equilibrium below cloud base:
471tdifexp=0.58
472tv0cin=7.0
473
474! facteurs lambda et mu dans Si=musi-lambda*T
475musi=1.0
476if (ok_nocinsat) then
477lambda_sursat = 0.0 ! no sursaturation effect
478endif           
479
480
481! diffusion dans le sol
482Kd=2.5e-9 ! m2/s   
483
484! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
485rh_cste_surf_cond=0.6
486T_cste_surf_cond=288.0
487
488!--------------------------------------------------------------
489! Parameters that depend on the nature of water isotopes:
490!--------------------------------------------------------------
491! ** constantes locales
492fac_enrichoce18=0.0005
493! on a alors tcor018=1+fac_enrichoce18
494! tcorD=1+fac_enrichoce18*8
495! tcorO17=1+fac_enrichoce18*0.528
496alpha_liq_sol_O18=1.00291 ! valeur de Lehmann & Siegenthaler, 1991,
497  ! Journal of Glaciology, vol 37, p 23
498talph1_O18=1137.
499talph2_O18=-0.4156
500talph3_O18=-2.0667E-3
501talps1_O18=11.839
502talps2_O18=-0.028244
503tkcin0_O18 = 0.006
504tkcin1_O18 = 0.000285
505tkcin2_O18 = 0.00082
506tdifrel_O18= 1./0.9723
507
508! rapport des ln(alphaeq) entre O18 et O17
509fac_coeff_eq17_liq=0.529 ! donné par Amaelle
510!      fac_coeff_eq17_ice=0.528 ! slope MWL
511fac_coeff_eq17_ice=0.529
512
513
514write(*,*) 'iso_O18,iso_HDO,iso_eau=',iso_O18,iso_HDO,iso_eau
515do 999 ixt = 1, niso
516write(*,*) 'iso_init 80: ixt=',ixt
517
518
519! -- kinetic factor for surface evaporation:
520! (cf: kcin = tkcin0                  if |V|<tv0cin
521!      kcin = tkcin1*|Vsurf| + tkcin2 if |V|>tv0cin )
522! (Rq: formula discontinuous for |V|=tv0cin... )       
523
524! -- main:
525if (ixt.eq.iso_HTO) then ! Tritium
526  tkcin0(ixt) = 0.01056
527  tkcin1(ixt) = 0.0005016
528  tkcin2(ixt) = 0.0014432
529  tnat(ixt)=0.
530  !toce(ixt)=2.2222E-8 ! corrigé par Alex Cauquoin
531  !toce(ixt)=1.0E-18 ! rapport 3H/1H ocean
532  toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
533  tcorr(ixt)=1.
534  tdifrel(ixt)=1./0.968
535  talph1(ixt)=46480.
536  talph2(ixt)=-103.87
537  talph3(ixt)=0.
538  talps1(ixt)=46480.
539  talps2(ixt)=-103.87
540  alpha_liq_sol(ixt)=1.
541  Rdefault(ixt)=0.0
542  Rmethox(ixt)=0.0
543  striso(ixt)='HTO'
544endif     
545if (ixt.eq.iso_O17) then ! Deuterium
546  pente_MWL=0.528
547!          tdifrel(ixt)=1./0.985452 ! donné par Amaelle
548  tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG
549!          fac_kcin=0.5145 ! donné par Amaelle
550  fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
551  tkcin0(ixt) = tkcin0_O18*fac_kcin
552  tkcin1(ixt) = tkcin1_O18*fac_kcin
553  tkcin2(ixt) = tkcin2_O18*fac_kcin
554  tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
555  toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
556  tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle           
557  talph1(ixt)=talph1_O18
558  talph2(ixt)=talph2_O18
559  talph3(ixt)=talph3_O18
560  talps1(ixt)=talps1_O18
561  talps2(ixt)=talps2_O18     
562  alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq
563  if (Rdefault_smow) then   
564        Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
565  else
566        Rdefault(ixt)=0.0
567  endif
568  Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
569  striso(ixt)='O17'
570endif
571
572if (ixt.eq.iso_O18) then ! Oxygene18
573  tkcin0(ixt) = tkcin0_O18
574  tkcin1(ixt) = tkcin1_O18
575  tkcin2(ixt) = tkcin2_O18
576  tnat(ixt)=2005.2E-6
577  toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
578  tcorr(ixt)=1.0+fac_enrichoce18
579  tdifrel(ixt)=tdifrel_O18
580  talph1(ixt)=talph1_O18
581  talph2(ixt)=talph2_O18
582  talph3(ixt)=talph3_O18
583  talps1(ixt)=talps1_O18
584  talps2(ixt)=talps2_O18
585  alpha_liq_sol(ixt)=alpha_liq_sol_O18   
586  if (Rdefault_smow) then   
587        Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
588  else
589        Rdefault(ixt)=0.0
590  endif
591  Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006   
592!       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
593  striso(ixt)='O18'
594  write(*,*) 'isotopes_mod 519: ixt,striso(ixt)=',ixt,striso(ixt)
595endif
596
597if (ixt.eq.iso_HDO) then ! Deuterium
598  pente_MWL=8.0
599!          fac_kcin=0.88
600  tdifrel(ixt)=1./0.9755
601  fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
602  tkcin0(ixt) = tkcin0_O18*fac_kcin
603  tkcin1(ixt) = tkcin1_O18*fac_kcin
604  tkcin2(ixt) = tkcin2_O18*fac_kcin
605  tnat(ixt)=155.76E-6
606  toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
607  tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL         
608  talph1(ixt)=24844.
609  talph2(ixt)=-76.248
610  talph3(ixt)=52.612E-3
611  talps1(ixt)=16288.
612  talps2(ixt)=-0.0934
613  !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
614  alpha_liq_sol(ixt)=1.0212
615  ! valeur de Lehmann & Siegenthaler, 1991, Journal of
616  ! Glaciology, vol 37, p 23
617  if (Rdefault_smow) then   
618    Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
619  else
620    Rdefault(ixt)=0.0
621  endif
622Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
623striso(ixt)='HDO'
624  write(*,*) 'isotopes_mod 548: ixt,striso(ixt)=',ixt,striso(ixt)
625endif
626
627!       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
628if (ixt.eq.iso_eau) then ! Oxygene16
629  tkcin0(ixt) = 0.0
630  tkcin1(ixt) = 0.0
631  tkcin2(ixt) = 0.0
632  tnat(ixt)=1.
633  toce(ixt)=tnat(ixt)
634  tcorr(ixt)=1.0
635  tdifrel(ixt)=1.
636  talph1(ixt)=0.
637  talph2(ixt)=0.
638  talph3(ixt)=0.
639  talps1(ixt)=0.
640  talph3(ixt)=0.
641  alpha_liq_sol(ixt)=1.
642  if (Rdefault_smow) then
643        Rdefault(ixt)=tnat(ixt)*1.0
644  else
645        Rdefault(ixt)=1.0
646  endif
647  Rmethox(ixt)=1.0
648  striso(ixt)='eau'
649endif
650
651999   continue
652
653! test de sensibilité:
654if (ok_nocinsfc) then ! no kinetic effect in sfc evaporation
655 do ixt=1,niso
656  tkcin0(ixt) = 0.0
657  tkcin1(ixt) = 0.0
658  tkcin2(ixt) = 0.0
659 enddo
660endif
661
662! fermeture fichier de paramètres
663close(unit=32)
664
665! nom des isotopes
666
667! verif
668write(*,*) 'iso_init 285: verif initialisation:'
669
670do ixt=1,niso
671  write(*,*) '* striso(',ixt,')=<'//striso(ixt)//'>'
672  write(*,*) 'tnat(',ixt,')=',tnat(ixt)
673!          write(*,*) 'alpha_liq_sol(',ixt,')=',alpha_liq_sol(ixt)
674!          write(*,*) 'tkcin0(',ixt,')=',tkcin0(ixt)
675!          write(*,*) 'tdifrel(',ixt,')=',tdifrel(ixt)
676enddo
677write(*,*) 'iso_init 69: lambda=',lambda_sursat
678write(*,*) 'iso_init 69: thumxt1=',thumxt1
679write(*,*) 'iso_init 69: h_land_ice=',h_land_ice
680write(*,*) 'iso_init 69: P_veg=',P_veg   
681           
682
683END SUBROUTINE iso_init
684
685! functions basiques mises ici pour éviter dépendances circulaires
686        !***********
687        function double_to_real(a)
688        implicit none
689        double precision a
690        real double_to_real
691
692        double_to_real=a
693
694        return
695        end function double_to_real
696
697        !***********
698        function real_to_double(a)
699        implicit none
700        real a
701        double precision real_to_double
702
703        real_to_double=a
704
705        return
706        end function real_to_double
707
708END MODULE isotopes_mod
709#endif
710
711
Note: See TracBrowser for help on using the repository browser.