[3466] | 1 | MODULE molvis_mod |
---|
| 2 | |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | |
---|
| 5 | CONTAINS |
---|
| 6 | |
---|
| 7 | SUBROUTINE molvis(ngrid,nlayer,ptimestep, |
---|
[1047] | 8 | & pplay,pplev,pt,pdteuv,pdtconduc |
---|
[690] | 9 | $ ,pvel,tsurf,zzlev,zzlay,zdvelmolvis) |
---|
[1047] | 10 | |
---|
| 11 | use conc_mod, only: cpnew, Akknew, rnew |
---|
[38] | 12 | IMPLICIT NONE |
---|
| 13 | |
---|
| 14 | c======================================================================= |
---|
| 15 | c |
---|
| 16 | c Molecular Viscosity Diffusion |
---|
| 17 | c |
---|
| 18 | c Based on conduction.F by N. Descamp, F. Forget 05/1999 |
---|
| 19 | c |
---|
| 20 | c modified by M. Angelats i Coll |
---|
| 21 | c |
---|
| 22 | c======================================================================= |
---|
| 23 | |
---|
| 24 | c----------------------------------------------------------------------- |
---|
| 25 | c declarations: |
---|
| 26 | c----------------------------------------------------------------------- |
---|
| 27 | |
---|
| 28 | c arguments: |
---|
| 29 | c ---------- |
---|
| 30 | |
---|
[1047] | 31 | integer,intent(in) :: ngrid ! number of atmospheric columns |
---|
| 32 | integer,intent(in) :: nlayer ! number of atmospheric layers |
---|
[38] | 33 | REAL ptimestep |
---|
[1047] | 34 | REAL pplay(ngrid,nlayer) |
---|
| 35 | REAL pplev(ngrid,nlayer+1) |
---|
| 36 | REAL zzlay(ngrid,nlayer) |
---|
| 37 | REAL zzlev(ngrid,nlayer+1) |
---|
| 38 | real pt(ngrid,nlayer) |
---|
| 39 | real tsurf(ngrid) |
---|
| 40 | REAL pvel(ngrid,nlayer) |
---|
| 41 | REAL pdvel(ngrid,nlayer) |
---|
| 42 | real pdteuv(ngrid,nlayer) |
---|
| 43 | real pdtconduc(ngrid,nlayer) |
---|
[38] | 44 | |
---|
[1047] | 45 | real zdvelmolvis(ngrid,nlayer) |
---|
[38] | 46 | |
---|
| 47 | c local: |
---|
| 48 | c ------ |
---|
| 49 | |
---|
[1047] | 50 | INTEGER l,ig, nz |
---|
| 51 | real Akk,phitop,fac, m, tmean |
---|
| 52 | REAL zvel(nlayer) |
---|
| 53 | real zt(nlayer) |
---|
| 54 | REAL alpha(nlayer) |
---|
| 55 | REAL lambda(nlayer) |
---|
| 56 | real muvol(nlayer) |
---|
| 57 | REAL C(nlayer) |
---|
| 58 | real D(nlayer) |
---|
| 59 | real den(nlayer) |
---|
| 60 | REAL pdvelm(nlayer) |
---|
| 61 | REAL zlay(nlayer) |
---|
| 62 | real zlev(nlayer+1) |
---|
[38] | 63 | |
---|
| 64 | c constants used locally |
---|
| 65 | c --------------------- |
---|
| 66 | c The atmospheric conductivity is a function of temperature T : |
---|
| 67 | c conductivity = Akk* T**skk |
---|
| 68 | c Molecular viscosity is related to thermal conductivity by: |
---|
| 69 | c conduc = 0.25*(9*gamma - 5)* Cv * molvis |
---|
| 70 | c where gamma = Cp/Cv. For dry air. |
---|
| 71 | |
---|
| 72 | |
---|
[1047] | 73 | REAL,PARAMETER :: skk=0.69 |
---|
[38] | 74 | |
---|
[1047] | 75 | REAL,PARAMETER :: velsurf =0.0 |
---|
[38] | 76 | |
---|
[1047] | 77 | logical,save :: firstcall=.true. |
---|
| 78 | |
---|
[2615] | 79 | !$OMP THREADPRIVATE(firstcall) |
---|
| 80 | |
---|
[38] | 81 | c----------------------------------------------------------------------- |
---|
| 82 | c calcul des coefficients alpha et lambda |
---|
| 83 | c----------------------------------------------------------------------- |
---|
| 84 | |
---|
| 85 | IF (firstcall) THEN |
---|
| 86 | ! write(*,*)'molvis: coeff of molecular viscosity Akk,skk,factor' |
---|
| 87 | ! write(*,*) Akk,skk,fac |
---|
| 88 | ! NB: Akk and fac are undefined at firstcall |
---|
| 89 | write(*,*)'molvis: coeff of molecular viscosity skk ', skk |
---|
| 90 | |
---|
| 91 | firstcall = .false. |
---|
| 92 | END IF |
---|
| 93 | |
---|
| 94 | ! Initialize phitop |
---|
| 95 | phitop=0.0 |
---|
| 96 | |
---|
[1047] | 97 | nz=nlayer |
---|
[38] | 98 | |
---|
| 99 | do ig=1,ngrid |
---|
| 100 | |
---|
| 101 | zt(1)=pt(ig,1)+(pdteuv(ig,1)+pdtconduc(ig,1))*ptimestep |
---|
| 102 | zvel(1)=pvel(ig,1) |
---|
| 103 | zlay(1)=zzlay(ig,1) |
---|
| 104 | zlev(1)=zzlev(ig,1) |
---|
| 105 | |
---|
| 106 | do l=2,nz |
---|
| 107 | zt(l)=pt(ig,l)+(pdteuv(ig,l)+pdtconduc(ig,l))*ptimestep |
---|
| 108 | zvel(l)=pvel(ig,l) |
---|
[3158] | 109 | zlay(l)=zzlay(ig,l) |
---|
| 110 | zlev(l)=zzlev(ig,l) |
---|
[38] | 111 | enddo |
---|
| 112 | |
---|
| 113 | zlev(nz+1)= zlev(nz)+10000. |
---|
| 114 | |
---|
| 115 | fac=0.25*(9.*cpnew(ig,1)-5.*(cpnew(ig,1)-rnew(ig,1))) |
---|
| 116 | Akk=Akknew(ig,1) |
---|
| 117 | lambda(1)=Akk*tsurf(ig)**skk/zlay(1)/fac |
---|
| 118 | c write(*,*) 'rnew(ig,nz) ',ig , rnew(ig,nz) |
---|
| 119 | |
---|
| 120 | DO l=2,nz |
---|
| 121 | fac=(9.*cpnew(ig,l)-5.*(cpnew(ig,l)-rnew(ig,l)))/4. |
---|
| 122 | Akk=Akknew(ig,l) |
---|
| 123 | lambda(l)=Akk/fac*zt(l)**skk/(zlay(l)-zlay(l-1)) |
---|
| 124 | ENDDO |
---|
| 125 | |
---|
| 126 | DO l=1,nz-1 |
---|
| 127 | muvol(l)=pplay(ig,l)/(rnew(ig,l)*zt(l)) |
---|
| 128 | alpha(l)=(muvol(l)/ptimestep)*(zlev(l+1)-zlev(l)) |
---|
| 129 | ENDDO |
---|
| 130 | muvol(nz)=pplay(ig,nz)/(rnew(ig,nz)*zt(nz)) |
---|
| 131 | alpha(nz)=(muvol(nz)/ptimestep)*(zlev(nz+1)-zlev(nz)) |
---|
| 132 | |
---|
| 133 | c-------------------------------------------------------------------- |
---|
| 134 | c |
---|
| 135 | c calcul des coefficients C et D |
---|
| 136 | c |
---|
| 137 | c------------------------------------------------------------------- |
---|
| 138 | |
---|
| 139 | den(1)=alpha(1)+lambda(2)+lambda(1) |
---|
| 140 | C(1)=lambda(1)*(velsurf-zvel(1))+lambda(2)*(zvel(2)-zvel(1)) |
---|
| 141 | C(1)=C(1)/den(1) |
---|
| 142 | D(1)=lambda(2)/den(1) |
---|
| 143 | |
---|
| 144 | DO l = 2,nz-1 |
---|
| 145 | den(l)=alpha(l)+lambda(l+1) |
---|
| 146 | den(l)=den(l)+lambda(l)*(1-D(l-1)) |
---|
| 147 | |
---|
| 148 | C(l) =lambda(l+1)*(zvel(l+1)-zvel(l)) |
---|
| 149 | $ +lambda(l)*(zvel(l-1)-zvel(l)+C(l-1)) |
---|
| 150 | C(l) =C(l)/den(l) |
---|
| 151 | |
---|
| 152 | D(l) =lambda(l+1) / den(l) |
---|
| 153 | ENDDO |
---|
| 154 | |
---|
| 155 | den(nz)=alpha(nz) + lambda(nz) * (1-D(nz-1)) |
---|
| 156 | C(nz)=C(nz-1)+zvel(nz-1)-zvel(nz) |
---|
| 157 | C(nz)=(C(nz)*lambda(nz)+phitop) / den(nz) |
---|
| 158 | |
---|
| 159 | c---------------------------------------------------------------------- |
---|
| 160 | c |
---|
| 161 | c calcul de la nouvelle pdvelm |
---|
| 162 | c |
---|
| 163 | c---------------------------------------------------------------------- |
---|
| 164 | |
---|
| 165 | DO l=1,nz |
---|
| 166 | pdvelm(l)=0. |
---|
| 167 | ENDDO |
---|
| 168 | pdvelm(nz)=C(nz) |
---|
| 169 | DO l=nz-1,1,-1 |
---|
| 170 | pdvelm(l)=C(l)+D(l)*pdvelm(l+1) |
---|
| 171 | ENDDO |
---|
| 172 | c----------------------------------------------------------------------- |
---|
| 173 | c |
---|
| 174 | c calcul de la tendance zdvelmolvis |
---|
| 175 | c |
---|
| 176 | c----------------------------------------------------------------------- |
---|
| 177 | |
---|
| 178 | DO l=1,nz |
---|
| 179 | zdvelmolvis(ig,l)=pdvelm(l)/ptimestep |
---|
| 180 | ENDDO |
---|
| 181 | |
---|
| 182 | ENDDO ! boucle sur ngrid |
---|
| 183 | |
---|
[3466] | 184 | END SUBROUTINE molvis |
---|
| 185 | |
---|
| 186 | END MODULE molvis_mod |
---|