| 1 | subroutine kastprof_fn(tsurf,temp,popsk,play,plev,qvap) |
|---|
| 2 | implicit none |
|---|
| 3 | |
|---|
| 4 | ! creates a kasting temperature profile |
|---|
| 5 | |
|---|
| 6 | #include "dimensions.h" |
|---|
| 7 | #include "dimphys.h" |
|---|
| 8 | #include "comcstfi.h" |
|---|
| 9 | #include "callkeys.h" |
|---|
| 10 | |
|---|
| 11 | REAL play(nlayermx) ! Pressure at the middle of the layers (Pa) |
|---|
| 12 | REAL plev(nlayermx+1) ! intermediate pressure levels (pa) |
|---|
| 13 | REAL psurf,tsurf |
|---|
| 14 | REAL tcond(nlayermx) |
|---|
| 15 | REAL qvap(nlayermx) |
|---|
| 16 | REAL temp(nlayermx) ! temperature at the middle of the layers |
|---|
| 17 | |
|---|
| 18 | real, parameter :: Rstar = 8.314 |
|---|
| 19 | real, parameter :: ptriple = 518000.0 |
|---|
| 20 | real, parameter :: latcond = 5.9e5 ! as in condense_co2 |
|---|
| 21 | real, parameter :: mn = 0.02801 ! molar mass of background gas (N2) |
|---|
| 22 | real, parameter :: mv = 0.04401 ! molar mass of background gas (CO2) |
|---|
| 23 | |
|---|
| 24 | real popsk(nlayermx) ! this is dodgy as in rest of code, zpopsk(ngridmx,nlayermx) |
|---|
| 25 | ! but kastprof requires 1D, so its ok I guess |
|---|
| 26 | |
|---|
| 27 | real ppco2 |
|---|
| 28 | real vmr_local,q_local,gamma,dqdlnT,dlnqdlnT,dlnpdlnT,dpdT,dln_p |
|---|
| 29 | |
|---|
| 30 | integer l,iter |
|---|
| 31 | |
|---|
| 32 | logical co2mixedcond |
|---|
| 33 | |
|---|
| 34 | co2mixedcond=.true. |
|---|
| 35 | |
|---|
| 36 | if(ngridmx.gt.1)then |
|---|
| 37 | print*,'DO NOT RUN KASTPROF IN 3D!!!' |
|---|
| 38 | call abort |
|---|
| 39 | end if |
|---|
| 40 | |
|---|
| 41 | ! 1st, we make pot. temperature constant everywhere |
|---|
| 42 | do l=1,nlayermx |
|---|
| 43 | temp(l)=TsurfK*popsk(l) |
|---|
| 44 | enddo |
|---|
| 45 | |
|---|
| 46 | ! then do all the other stuff |
|---|
| 47 | vmr_local = 1.-n2mixratio |
|---|
| 48 | q_local = vmr_local*mv/mn |
|---|
| 49 | |
|---|
| 50 | qvap(1) = q_local |
|---|
| 51 | |
|---|
| 52 | ppco2=vmr_local*play(1) |
|---|
| 53 | |
|---|
| 54 | if(play(1).lt.ptriple)then ! |
|---|
| 55 | tcond(1) = (-3167.8)/(log(.01*ppco2)-23.23) ! Fanale's formula |
|---|
| 56 | else |
|---|
| 57 | tcond(1) = 684.2-92.3*log(ppco2)+4.32*log(ppco2)**2 ! liquid-vapour (CRC handbook 2003 data) |
|---|
| 58 | endif |
|---|
| 59 | |
|---|
| 60 | do l=1,nlayermx-1 |
|---|
| 61 | |
|---|
| 62 | ppco2=vmr_local*play(l) ! this shld be implicit i.e. vmr_local(l+1) ?? |
|---|
| 63 | |
|---|
| 64 | if(play(l).le.ptriple)then ! ERROR |
|---|
| 65 | tcond(l) = (-3167.8)/(log(.01*ppco2)-23.23) ! Fanale's formula |
|---|
| 66 | else |
|---|
| 67 | tcond(l) = 684.2-92.3*log(ppco2)+4.32*log(ppco2)**2 ! liquid-vapour (CRC handbook 2003 data) |
|---|
| 68 | endif |
|---|
| 69 | |
|---|
| 70 | if(co2mixedcond)then |
|---|
| 71 | if (temp(l) .le. Tstrat)then |
|---|
| 72 | temp(l) = Tstrat ! isothermal stratosphere |
|---|
| 73 | temp(l+1) = Tstrat |
|---|
| 74 | qvap(l+1) = qvap(l) |
|---|
| 75 | !print*,'strato' |
|---|
| 76 | |
|---|
| 77 | elseif(temp(l) .le. tcond(l))then |
|---|
| 78 | |
|---|
| 79 | dln_p = log(play(l+1))-log(play(l)) |
|---|
| 80 | gamma = cpp - q_local*latcond/temp(l) |
|---|
| 81 | dqdlnT = (mv*latcond/(mn*temp(l)) - gamma) & |
|---|
| 82 | / (latcond/temp(l) + Rstar/(mn*q_local)) |
|---|
| 83 | dlnqdlnT = (1./q_local)*dqdlnT |
|---|
| 84 | dlnpdlnT = mv*latcond/(Rstar*temp(l)) - dlnqdlnT/(1. + q_local*mn/mv) |
|---|
| 85 | dpdT = (play(l)/temp(l))*dlnpdlnT |
|---|
| 86 | q_local = q_local + dqdlnT*dln_p/dlnpdlnT |
|---|
| 87 | |
|---|
| 88 | vmr_local = q_local*mn/mv |
|---|
| 89 | |
|---|
| 90 | |
|---|
| 91 | if(1.eq.1)then |
|---|
| 92 | print*,'l =',l |
|---|
| 93 | print*,'gamma =',gamma |
|---|
| 94 | print*,'gamma bit =',q_local*latcond/temp(l) |
|---|
| 95 | print*,'cpp =',cpp |
|---|
| 96 | print*,'latcond =',latcond |
|---|
| 97 | print*,'q_local =',q_local |
|---|
| 98 | print*,'T =',temp(l) |
|---|
| 99 | print*,'p =',play(l) |
|---|
| 100 | print*,'dqdlnT =',dqdlnT |
|---|
| 101 | print*,'dlnqdlnT =',dlnqdlnT |
|---|
| 102 | print*,'dlnpdlnT =',dlnpdlnT |
|---|
| 103 | print*,'dpdT =',dpdT |
|---|
| 104 | call abort |
|---|
| 105 | endif |
|---|
| 106 | |
|---|
| 107 | qvap(l+1) = q_local |
|---|
| 108 | temp(l+1) = temp(l) + temp(l)*dln_p/(dlnpdlnT) |
|---|
| 109 | !print*,'CO2co' |
|---|
| 110 | else |
|---|
| 111 | !print*,'tropo' |
|---|
| 112 | |
|---|
| 113 | |
|---|
| 114 | ! CO2 is not condensing; therefore dry convection |
|---|
| 115 | qvap(l+1) = qvap(l) |
|---|
| 116 | endif |
|---|
| 117 | |
|---|
| 118 | else |
|---|
| 119 | |
|---|
| 120 | if (temp(l) .le. Tstrat)then |
|---|
| 121 | temp(l) = Tstrat |
|---|
| 122 | ! isothermal stratosphere |
|---|
| 123 | elseif(temp(l) .le. tcond(l))then |
|---|
| 124 | temp(l)=tcond(l) |
|---|
| 125 | ! pure CO2 condensation |
|---|
| 126 | else |
|---|
| 127 | ! CO2 is not condensing; therefore dry convection |
|---|
| 128 | endif |
|---|
| 129 | |
|---|
| 130 | endif |
|---|
| 131 | |
|---|
| 132 | enddo ! vertical loop |
|---|
| 133 | |
|---|
| 134 | temp(nlayermx) = Tstrat |
|---|
| 135 | qvap(nlayermx) = qvap(nlayermx-1) |
|---|
| 136 | |
|---|
| 137 | print*,'---fin---' |
|---|
| 138 | |
|---|
| 139 | call abort |
|---|
| 140 | |
|---|
| 141 | return |
|---|
| 142 | end subroutine kastprof_fn |
|---|
