Discussion :

[leovenus]

Bonjour!

J'ai écrit une petite partie de code me permettant de modifier la
température d'une couche de neige dans un modèle 1D RCM.J'ai écrit en parallèle le même code en fortran et en matlab (pour faire
mes vérifications). Matlab marche, Fortran ne marche pas... Dans les codes
ci-dessous, T0 a la bonne valeur dans les deux cas, mais quand je veux
assigner T0 à at(19), il me note que at(19)=0.00000. Est-ce un problème de syntaxe?
(J'ai tout vérifié, toutes les valeurs sont identiques sauf at(19)=T0)

Lucie

FORTRAN : (tout ce qui se trouve en-dessous de * Boucle while etc.)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
 
      subroutine tempchng(pa,p,cp,at,dg,df,g,step,s1,f,ch,ce,
     xkap,abt,c,at1)
      real pa(20),at(20),dg(20),df(20),f(20),g(20),ch(19)
      real p(20),cp(20),ce(20),c(20),at1(20)
      real dt,s1,dr,step,e1,e,q,l
      integer kap
      real corr,T0,K,Tp,fct,fctPrime,sigma,d
      abt=0.0
      et=0.0
      s1=0.0
      cp(19)=5*4.186e6/86400
* MODIFICATION : prise en compte de la glace plutôt que de l'eau
* Variation de at(19) selon cette modification
      corr=2.0
      T0=at(19)
      K=20.0
*     Température de glace profonde : -50°C
      Tp=223.15
      fct=0.0
      fctPrime=0.0
      sigma=5.67e-8
      d=step/cp(19)
*
      do 3400 j=1,18
*       if (j.eq.kap) then
*        h=1.0
*       end if
       e1=esat(at1(j)+1)
       e=esat(at1(j))
*     calculate modified heat capacity Manabe & Wetherald 1967?
      dr=rwat(pa(j),e1,j,kap)-rwat(pa(j),e,j,kap)
      l=2510-2.38*(at(j)-273)
       c(j)=.622*l*l*rwat(pa(j),e,j,kap)/(1.005*.287*at1(j)**2)
      cp(j)=(1.0+c(j))*1.038165e7*(p(j+1)-p(j))/86400
      dt=step*(dg(j)+((ce(j))*(1.0+c(j))/step)-df(j))
     x/(1.0+c(j))
       et=et+cp(j)*(ce(j))
       s1=s1+abs(dt+ch(j))
       at(j)=at(j)+dt
3400  continue
       et=et/step
        ce(19)=-et*step/cp(19)
*
* Boucle while en fortran
666   if(corr.gt.0.1) then
       fct=T0+d*(sigma*T0**4+T0)-at(19)-d*(g(19)-et-f(19)+K*Tp)
       fctPrime=4*step*sigma*T0**3/cp(19)+1+d
       corr= -(fct/fctPrime)
       T0=T0+corr
       goto 666
      endif
      at(19)=T0
*      write(40,*),at(19)
       return
        end

Matlab :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
 
function at = tempchng()
pa= [2.2719479E-03,1.9675927E-02,5.2512009E-02,9.8722570E-02,0.1562500,...
    0.2230367,0.2970250,0.3761574,0.4583762,0.5416238,0.6238427,0.7029751,...
    0.7769634,0.8437501,0.9012776,0.9474881,0.9803241,0.9977280,1.0];
p=[0,9.0020578E-03,3.4465022E-02,7.4331276E-02,0.1265432,0.1890432,...
    0.2597737,0.3366770,0.4176955,0.5007716,0.5838478,0.6648663,0.7417696...
    0.8125000,0.8750000,0.9272119,0.9670781,0.9925410,1.0];
cp=zeros(19,1);
at=[259.1344,236.0187,209.7579,207.8609,207.3105,221.0353,233.3348,243.9981,...
    253.3060,261.4459,268.5537,274.7156,279.9859,284.4055,287.9909,290.7393,...
    292.6271,293.6071,293.7351];
dg=[1.335051,0.5880523,0.2164474,0.1048099,6.9614537E-02,9.3966797E-02,0.1313506...    
    0.1406944,0.1491358,0.1625618,0.1020992,9.1658607E-02, 9.1098346E-02,...
    9.0158798E-02,8.9044750E-02,8.7956138E-02,8.7107100E-02,8.6886391E-02,0];
df=[4.777882,1.720254,0.5534142,0.1319819,0.1397153,1.002216,1.745454...
    1.935682,1.865595,1.826197,4.407408,0.5994710,0.7716145,0.8345697...
    0.8792982,0.9926330,1.144125,1.336380,0];
g=[42.52134,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,12.35275];
step=1;
s1=0;
f=[246.8946,0,0,0,0,-231.935,0,0,0,0,0,0,0,0,0,0,0,-75.35191,346.7395];
ch=zeros(19,1);
ce=zeros(29,1);
kap=11;
abt=0;
c=zeros(19,1);
at1=[259.1344,236.0187,209.7579,207.8609,207.3105,221.0353,233.3348,...
    243.9981,253.306,261.4459,268.5537,274.7156,279.9859,284.4055,287.9909,...
    290.7393,292.6271,293.6071,293.7351];
et=0;
corr=2;
T0=at(19);
K=20;
Tp=223.15;
fct=0;
fctPr=0;
cp(19)=5*4.186e6/86400;
 
for j=1:18
    if(j==kap)
        h=1;
    end
    e1=(6.11/1012.34)*exp((.622*(2510.-2.38*((at1(j)+1)-273))/.287)*((at1(j)+1)-273)/((at1(j)+1)*273));
    e=(6.11/1012.34)*exp((.622*(2510.-2.38*(at1(j)-273))/.287)*(at1(j)-273)/(at1(j)*273));
% calculate modified heat capacity Manabe & Wetherald 1967?
    dr=rwat(pa(j),e1,j,kap)-rwat(pa(j),e,j,kap);
    l=2510-2.38*(at(j)-273);
    c(j)=.622*l^2*rwat(pa(j),e,j,kap)/(1.005*.287*at1(j)^2);
    cp(j)=(1+c(j))*1.038165e7*(p(j+1)-p(j))/86400;
    dt=step*(dg(j)+((ce(j))*(1+c(j))/step)-df(j))/(1+c(j));
    et=et+cp(j)*ce(j);
    s1=s1+abs(dt+ch(j));
    at(j)=at(j)+dt;
end
sigma=5.67E-8;
et=et/step
ce(19)=-et*step/cp(19);
while abs(corr)>0.1
    fct=T0+(step/cp(19))*(sigma*T0^4+T0)-at(19)-(step/cp(19))*(g(19)-et-f(19)+K*Tp);
    fctPr=4*step*sigma*T0^3/cp(19)+1+step/cp(19);
    corr= -fct/fctPr;
    T0=T0+corr;
end
at(19)=T0;
end
 
%pour "rwat":
    function rwat= rwat(p,e,j,kap)
        if(j==kap)
            h=.77*(p-.02)/.98;
        else
            h=.77*(p-.02)/.98;
        end
        rwat=.622*h*e/(p-h*e);
        if (rwat<3e-6)
            rwat=3e-6;
        end
    end

[Sylvain Bergeron]

Je n'ai pas vraiment regardé le code, sauf que je remarque qu'il n'y a pas d'instruction « implicit none ». À mon avis, tenter de chercher une erreur comme celle que tu observes sans d'abord ajouter un « implicit none » est une perte de temps...