| 1 | subroutine N2broadprof_fn(tsurf,temp,popsk,play,plev,qvap) |
|---|
| 2 | implicit none |
|---|
| 3 | |
|---|
| 4 | ! creates a temperature profile for an N2-CO2 atmosphere |
|---|
| 5 | ! a fine resolution grid is used for accuracy |
|---|
| 6 | |
|---|
| 7 | #include "dimensions.h" |
|---|
| 8 | #include "dimphys.h" |
|---|
| 9 | #include "comcstfi.h" |
|---|
| 10 | #include "callkeys.h" |
|---|
| 11 | |
|---|
| 12 | REAL play(nlayermx) ! Pressure at the middle of the layers (Pa) |
|---|
| 13 | REAL plev(nlayermx+1) ! Intermediate pressure levels (pa) |
|---|
| 14 | REAL psurf,tsurf |
|---|
| 15 | ! REAL tcond(nlayermx) |
|---|
| 16 | REAL qvap(nlayermx) |
|---|
| 17 | REAL temp(nlayermx) ! temperature at the middle of the layers |
|---|
| 18 | |
|---|
| 19 | real, parameter :: Rstar = 8.314 |
|---|
| 20 | real, parameter :: ptriple = 518000.0 |
|---|
| 21 | real, parameter :: latcond = 5.9e5 ! as in condense_co2 |
|---|
| 22 | real, parameter :: mn = 0.02801 ! molar mass of background gas (N2) |
|---|
| 23 | real, parameter :: mv = 0.04401 ! molar mass of background gas (CO2) |
|---|
| 24 | |
|---|
| 25 | real popsk(nlayermx) ! this is dodgy as in rest of code, zpopsk(ngridmx,nlayermx) |
|---|
| 26 | ! but kastprof requires 1D, so its ok I guess |
|---|
| 27 | |
|---|
| 28 | real gamma,dqdlnT,dlnqdlnT,dlnpdlnT,dpdT,dln_p |
|---|
| 29 | |
|---|
| 30 | integer l,iter |
|---|
| 31 | |
|---|
| 32 | real sigma_loc,q_loc,vmr_loc,ppvar_loc,T_loc,p_loc |
|---|
| 33 | real Tcond |
|---|
| 34 | real Told,qold,pold |
|---|
| 35 | |
|---|
| 36 | integer, parameter :: Niter = 2000 |
|---|
| 37 | |
|---|
| 38 | if(ngridmx.gt.1)then |
|---|
| 39 | print*,'DO NOT RUN KASTPROF IN 3D!!!' |
|---|
| 40 | call abort |
|---|
| 41 | end if |
|---|
| 42 | |
|---|
| 43 | l=1 |
|---|
| 44 | |
|---|
| 45 | vmr_loc = 0.05 ! value at ground |
|---|
| 46 | q_loc = vmr_loc*mv/mn |
|---|
| 47 | T_loc = Tsurf |
|---|
| 48 | |
|---|
| 49 | dln_p = log(pceil/plev(1))/Niter ! please check |
|---|
| 50 | do iter=1,Niter |
|---|
| 51 | |
|---|
| 52 | sigma_loc = (pceil/plev(1))**(real(iter)/real(Niter)) |
|---|
| 53 | p_loc = plev(1)*sigma_loc |
|---|
| 54 | |
|---|
| 55 | ppvar_loc = vmr_loc*p_loc |
|---|
| 56 | if(ppvar_loc.lt.ptriple)then |
|---|
| 57 | Tcond = (-3167.8)/(log(.01*ppvar_loc)-23.23) ! Fanale's formula |
|---|
| 58 | else |
|---|
| 59 | Tcond = 684.2-92.3*log(ppvar_loc)+4.32*log(ppvar_loc)**2 ! liquid-vapour (CRC handbook 2003 data) |
|---|
| 60 | endif |
|---|
| 61 | |
|---|
| 62 | if (T_loc.le.Tstrat)then |
|---|
| 63 | T_loc = Tstrat ! isothermal stratosphere |
|---|
| 64 | elseif(T_loc.le.Tcond)then |
|---|
| 65 | |
|---|
| 66 | gamma = cpp - q_loc*latcond/T_loc |
|---|
| 67 | dqdlnT = (mv*latcond/(mn*T_loc) - gamma) & |
|---|
| 68 | / (latcond/T_loc + Rstar/(mn*q_loc)) |
|---|
| 69 | dlnqdlnT = (1./q_loc)*dqdlnT |
|---|
| 70 | dlnpdlnT = mv*latcond/(Rstar*T_loc) - dlnqdlnT/(1. + q_loc*mn/mv) |
|---|
| 71 | dpdT = (p_loc/T_loc)*dlnpdlnT |
|---|
| 72 | q_loc = q_loc + dqdlnT*dln_p/dlnpdlnT |
|---|
| 73 | vmr_loc = q_loc*mn/mv ! needed for Tcond calc. |
|---|
| 74 | T_loc = T_loc*(1 + dln_p/dlnpdlnT) |
|---|
| 75 | print*,T_loc |
|---|
| 76 | else |
|---|
| 77 | ! nothing doing; so dry convection |
|---|
| 78 | T_loc = tsurf*(p_loc/plev(1))**rcp |
|---|
| 79 | endif |
|---|
| 80 | |
|---|
| 81 | ! update true gridpoints when we reach them |
|---|
| 82 | if((p_loc.le.play(l)) .and. (l.le.nlayermx))then |
|---|
| 83 | temp(l)=T_loc |
|---|
| 84 | qvap(l)=q_loc |
|---|
| 85 | ! do an avg you lazy bastard |
|---|
| 86 | ! temp(l)=T_loc*(p_loc-play(l))/(p_loc-pold)+Told*(play(l)-pold)/(p_loc-pold) |
|---|
| 87 | ! qvap(l)=q_loc*(p_loc-play(l))/(p_loc-pold)+qold*(play(l)-pold)/(p_loc-pold) |
|---|
| 88 | |
|---|
| 89 | ! f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1) |
|---|
| 90 | |
|---|
| 91 | ! print*,'l=',l |
|---|
| 92 | ! print*,'temp(l)=',temp(l) |
|---|
| 93 | ! print*,'qvap(l)=',qvap(l) |
|---|
| 94 | |
|---|
| 95 | |
|---|
| 96 | qvap(l)=0.05*mv/mn ! suface amount everywhere to test |
|---|
| 97 | |
|---|
| 98 | l=l+1 |
|---|
| 99 | endif |
|---|
| 100 | |
|---|
| 101 | Told=T_loc |
|---|
| 102 | qold=q_loc |
|---|
| 103 | pold=p_loc |
|---|
| 104 | |
|---|
| 105 | enddo ! vertical loop |
|---|
| 106 | |
|---|
| 107 | print*,'---fin---' |
|---|
| 108 | |
|---|
| 109 | return |
|---|
| 110 | end subroutine N2broadprof_fn |
|---|
