source: trunk/LMDZ.MARS/libf/aeronomars/conduction.F

Last change on this file was 3466, checked in by emillour, 2 months ago

Mars PCM:
More tidying in aeronomars:

  • remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
  • turn aeronomars routines to modules, for those which aren't in modules yet.

EM

File size: 4.9 KB
RevLine 
[3466]1      MODULE conduction_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6     
7      SUBROUTINE conduction(ngrid,nlayer,ptimestep,pplay,pplev,pt,pdt,
[38]8     $                   tsurf,zzlev,zzlay,zdtconduc)
9
[1047]10      use conc_mod, only: Akknew, rnew, cpnew
[38]11      IMPLICIT NONE
12
13c=======================================================================
14c
15c   Molecular thermal conduction
16c   
17c   N. Descamp, F. Forget 05/1999
18c
19c=======================================================================
20
21c-----------------------------------------------------------------------
22c   declarations:
23c-----------------------------------------------------------------------
24
25c   arguments:
26c   ----------
27
[1047]28      integer,intent(in) :: ngrid ! number of atmospheric columns
29      integer,intent(in) :: nlayer ! number of atmospheric layers
30      real,intent(in) :: ptimestep
31      REAL,intent(in) :: pplay(ngrid,nlayer)
32      real,intent(in) :: pplev(ngrid,nlayer+1)
33      REAL,intent(in) :: zzlay(ngrid,nlayer)
34      real,intent(in) :: zzlev(ngrid,nlayer+1)
35      REAL,intent(in) :: pt(ngrid,nlayer)
36      real,intent(in) :: pdt(ngrid,nlayer)
37      real,intent(in) :: tsurf(ngrid)
[38]38
[1047]39      real,intent(out) :: zdtconduc(ngrid,nlayer)
[38]40
41c   local:
42c   ------
43
44      INTEGER i,ig,l
45      real Akk
46      real,save :: phitop
47      real m,tmean
[1047]48      REAL alpha(nlayer)
49      real zt(nlayer)
50      REAL lambda(nlayer)
51      real muvol(nlayer)
52      REAL C(nlayer)
53      real D(nlayer)
54      real den(nlayer)
55      REAL pdtc(nlayer)
56      real zlay(nlayer)
57      real zlev(nlayer+1)
[38]58
59c   constants used locally
60c    ---------------------
61c     The atmospheric conductivity is a function of temperature T :
62c      conductivity = Akk* T**skk
63      REAL,PARAMETER :: skk=0.69
64     
65      logical,save :: firstcall=.true.
66
[2615]67!$OMP THREADPRIVATE(phitop,firstcall)
68
[38]69c-----------------------------------------------------------------------
70c   calcul des coefficients alpha et lambda
71c-----------------------------------------------------------------------
72
73      IF (firstcall) THEN
74!        write (*,*)'conduction: coeff to compute molecular',
75!     &             ' conductivity Akk,skk'
76!        write(*,*) Akk,skk
77! NB: Akk is undefined at this stage
78        write (*,*)'conduction: coeff to compute molecular',
79     &             ' conductivity skk = ', skk
80
81! Initialize phitop
82        phitop=0.0
83       
84        firstcall = .false.
85      ENDIF ! of IF (firstcall)
86
87      do ig=1,ngrid
88
89        zt(1)=pt(ig,1)+pdt(ig,1)*ptimestep
90        zlay(1)=zzlay(ig,1)
91        zlev(1)=zzlev(ig,1)
92     
[1047]93        do i=2,nlayer
[38]94
95          zt(i)=pt(ig,i)+pdt(ig,i)*ptimestep
[3158]96          zlay(i)=zzlay(ig,i)
97          zlev(i)=zzlev(ig,i)
[38]98        enddo
[3158]99           
[1047]100        zlev(nlayer+1)= zlev(nlayer)+10000.
[38]101
102        Akk=Akknew(ig,1)
103        lambda(1) = Akk*tsurf(ig)**skk/zlay(1)   
104
[1047]105        DO i = 2 , nlayer
[38]106          Akk=Akknew(ig,i)
107          lambda(i)=Akk*zt(i)**skk/(zlay(i)-zlay(i-1))
108        ENDDO
[1047]109        DO i=1,nlayer-1
[38]110          muvol(i)=pplay(ig,i)/(rnew(ig,i)*zt(i))
111          alpha(i)=cpnew(ig,i)*(muvol(i)/ptimestep)
112     $                        *(zlev(i+1)-zlev(i))
113        ENDDO
114
[1047]115        muvol(nlayer)=pplay(ig,nlayer)/(rnew(ig,nlayer)*zt(nlayer))
116        alpha(nlayer)=cpnew(ig,i)*(muvol(nlayer)/ptimestep)
117     $                       *(zlev(nlayer+1)-zlev(nlayer))
[38]118
119c--------------------------------------------------------------------
120c
121c     calcul des coefficients C et D
122c
123c-------------------------------------------------------------------
124
125        den(1)=alpha(1)+lambda(2)+lambda(1)
126        C(1)=lambda(1)*(tsurf(ig)-zt(1))+lambda(2)*(zt(2)-zt(1))
127        C(1)=C(1)/den(1)             
128        D(1)=lambda(2)/den(1)           
129   
[1047]130        DO i = 2,nlayer-1
[38]131          den(i)=alpha(i)+lambda(i+1)
132          den(i)=den(i)+lambda(i)*(1-D(i-1))
133           
134          C(i) =lambda(i+1)*(zt(i+1)-zt(i))
135     $         +lambda(i)*(zt(i-1)-zt(i)+C(i-1))   
136          C(i) =C(i)/den(i)           
137
138          D(i) =lambda(i+1) / den(i)
139        ENDDO
140
[1047]141        den(nlayer)=alpha(nlayer) + lambda(nlayer) * (1-D(nlayer-1))
142        C(nlayer)=C(nlayer-1)+zt(nlayer-1)-zt(nlayer)
143        C(nlayer)=(C(nlayer)*lambda(nlayer)+phitop) / den(nlayer)
[38]144               
145c----------------------------------------------------------------------
146c
147c      calcul de la nouvelle temperature ptconduc
148c
149c----------------------------------------------------------------------
150
[1047]151        DO i=1,nlayer
[38]152          pdtc(i)=0.
153        ENDDO
[1047]154        pdtc(nlayer)=C(nlayer)
155        DO i=nlayer-1,1,-1
[38]156          pdtc(i)=C(i)+D(i)*pdtc(i+1)
157        ENDDO
158c-----------------------------------------------------------------------
159c
160c     calcul de la tendance zdtconduc
161c
162c-----------------------------------------------------------------------
163   
[1047]164        DO i=1,nlayer
[38]165          zdtconduc(ig,i)=pdtc(i)/ptimestep
166        ENDDO
167
168      enddo ! of do ig=1,ngrid
169
[3466]170      END SUBROUTINE conduction
171     
172      END MODULE conduction_mod
Note: See TracBrowser for help on using the repository browser.