| 1 | program kcm1d |
|---|
| 2 | |
|---|
| 3 | use radinc_h, only: NAERKIND |
|---|
| 4 | use radcommon_h, only: L_NSPECTI, L_NSPECTV, sigma |
|---|
| 5 | use watercommon_h, only: mH2O |
|---|
| 6 | use ioipsl_getincom |
|---|
| 7 | |
|---|
| 8 | implicit none |
|---|
| 9 | |
|---|
| 10 | !================================================================== |
|---|
| 11 | ! |
|---|
| 12 | ! Purpose |
|---|
| 13 | ! ------- |
|---|
| 14 | ! Run the universal model radiative transfer once in a 1D column. |
|---|
| 15 | ! Useful for climate evolution studies etc. |
|---|
| 16 | ! |
|---|
| 17 | ! It can be compiled with a command like (e.g. for 25 layers): |
|---|
| 18 | ! "makegcm -p std -d 25 kcm1d" |
|---|
| 19 | ! It requires the files "callphys.def", "gases.def" |
|---|
| 20 | ! "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def" |
|---|
| 21 | ! |
|---|
| 22 | ! Authors |
|---|
| 23 | ! ------- |
|---|
| 24 | ! R. Wordsworth |
|---|
| 25 | ! |
|---|
| 26 | !================================================================== |
|---|
| 27 | |
|---|
| 28 | #include "dimensions.h" |
|---|
| 29 | #include "dimphys.h" |
|---|
| 30 | #include "callkeys.h" |
|---|
| 31 | #include "comcstfi.h" |
|---|
| 32 | #include "planete.h" |
|---|
| 33 | #include "comsaison.h" |
|---|
| 34 | #include "control.h" |
|---|
| 35 | #include "advtrac.h" |
|---|
| 36 | |
|---|
| 37 | ! -------------------------------------------------------------- |
|---|
| 38 | ! Declarations |
|---|
| 39 | ! -------------------------------------------------------------- |
|---|
| 40 | |
|---|
| 41 | integer nlayer,nlevel,nq |
|---|
| 42 | integer ilay,ilev,iq,iw,iter |
|---|
| 43 | real play(nlayermx) ! Pressure at the middle of the layers [Pa] |
|---|
| 44 | real zlay(nlayermx) ! Altitude at middle of the layers [km] |
|---|
| 45 | real plev(nlayermx+1) ! Intermediate pressure levels [Pa] |
|---|
| 46 | real temp(nlayermx) ! temperature at the middle of the layers [K] |
|---|
| 47 | real q(nlayermx,nqmx) ! tracer mixing ratio [kg/kg] |
|---|
| 48 | real vmr(nlayermx,nqmx) ! tracer mixing ratio [mol/mol] |
|---|
| 49 | real qsurf(nqmx) ! tracer surface budget [kg/kg] ???? |
|---|
| 50 | real psurf,psurf_n,tsurf |
|---|
| 51 | |
|---|
| 52 | real emis, albedo |
|---|
| 53 | |
|---|
| 54 | real muvar(nlayermx+1) |
|---|
| 55 | |
|---|
| 56 | real dtsw(nlayermx) ! heating rate (K/s) due to SW |
|---|
| 57 | real dtlw(nlayermx) ! heating rate (K/s) due to LW |
|---|
| 58 | real fluxsurf_lw ! incident LW flux to surf (W/m2) |
|---|
| 59 | real fluxtop_lw ! outgoing LW flux to space (W/m2) |
|---|
| 60 | real fluxsurf_sw ! incident SW flux to surf (W/m2) |
|---|
| 61 | real fluxabs_sw ! SW flux absorbed by planet (W/m2) |
|---|
| 62 | real fluxtop_dn ! incident top of atmosphere SW flux (W/m2) |
|---|
| 63 | |
|---|
| 64 | ! not used |
|---|
| 65 | real reffrad(nlayermx,naerkind) |
|---|
| 66 | real nueffrad(nlayermx,naerkind) |
|---|
| 67 | real cloudfrac(nlayermx) |
|---|
| 68 | real totcloudfrac |
|---|
| 69 | real tau_col |
|---|
| 70 | |
|---|
| 71 | real dTstrat |
|---|
| 72 | real aerosol(nlayermx,naerkind) ! aerosol tau (kg/kg) |
|---|
| 73 | real OLR_nu(L_NSPECTI) |
|---|
| 74 | real GSR_nu(L_NSPECTV) |
|---|
| 75 | real Eatmtot |
|---|
| 76 | |
|---|
| 77 | integer ierr |
|---|
| 78 | logical firstcall,lastcall |
|---|
| 79 | |
|---|
| 80 | ! -------------- |
|---|
| 81 | ! Initialisation |
|---|
| 82 | ! -------------- |
|---|
| 83 | |
|---|
| 84 | pi=2.E+0*asin(1.E+0) |
|---|
| 85 | |
|---|
| 86 | reffrad(:,:) = 0.0 |
|---|
| 87 | nueffrad(:,:) = 0.0 |
|---|
| 88 | cloudfrac(:) = 0.0 |
|---|
| 89 | totcloudfrac = 0.0 |
|---|
| 90 | |
|---|
| 91 | ! Load parameters from "run.def" |
|---|
| 92 | ! ------------------------------- |
|---|
| 93 | |
|---|
| 94 | ! check if 'run1d.def' file is around (otherwise reading parameters |
|---|
| 95 | ! from callphys.def via getin() routine won't work.) |
|---|
| 96 | open(90,file='run.def',status='old',form='formatted',& |
|---|
| 97 | iostat=ierr) |
|---|
| 98 | if (ierr.ne.0) then |
|---|
| 99 | write(*,*) 'Cannot find required file "run.def"' |
|---|
| 100 | write(*,*) ' (which should contain some input parameters' |
|---|
| 101 | write(*,*) ' along with the following line:' |
|---|
| 102 | write(*,*) ' INCLUDEDEF=callphys.def' |
|---|
| 103 | write(*,*) ' )' |
|---|
| 104 | write(*,*) ' ... might as well stop here ...' |
|---|
| 105 | stop |
|---|
| 106 | else |
|---|
| 107 | close(90) |
|---|
| 108 | endif |
|---|
| 109 | |
|---|
| 110 | nlayer=nlayermx |
|---|
| 111 | nlevel=nlayer+1 |
|---|
| 112 | |
|---|
| 113 | psurf_n=100000. ! default value for psurf |
|---|
| 114 | print*,'Dry surface pressure (Pa)?' |
|---|
| 115 | call getin("psurf",psurf_n) |
|---|
| 116 | write(*,*) " psurf = ",psurf_n |
|---|
| 117 | |
|---|
| 118 | tsurf=300.0 ! default value for tsurf |
|---|
| 119 | print*,'Surface temperature (K)?' |
|---|
| 120 | call getin("tref",tsurf) |
|---|
| 121 | write(*,*) " tsurf = ",tsurf |
|---|
| 122 | |
|---|
| 123 | g=10.0 ! default value for g |
|---|
| 124 | print*,'Gravity ?' |
|---|
| 125 | call getin("g",g) |
|---|
| 126 | write(*,*) " g = ",g |
|---|
| 127 | |
|---|
| 128 | periastr = 1.0 |
|---|
| 129 | apoastr = 1.0 |
|---|
| 130 | print*,'Periastron (AU)?' |
|---|
| 131 | call getin("periastr",periastr) |
|---|
| 132 | write(*,*) "periastron = ",periastr |
|---|
| 133 | dist_star = periastr |
|---|
| 134 | fract = 0.5 |
|---|
| 135 | print*,'Apoastron (AU)?' |
|---|
| 136 | call getin("apoastr",apoastr) |
|---|
| 137 | write(*,*) "apoastron = ",apoastr |
|---|
| 138 | |
|---|
| 139 | periastr = periastr*149.598 ! AU to Mkm |
|---|
| 140 | apoastr = apoastr*149.598 ! AU to Mk |
|---|
| 141 | |
|---|
| 142 | albedo=0.2 ! default value for albedo |
|---|
| 143 | print*,'Albedo of bare ground?' |
|---|
| 144 | call getin("albedo",albedo) |
|---|
| 145 | write(*,*) " albedo = ",albedo |
|---|
| 146 | |
|---|
| 147 | emis=1.0 ! default value for emissivity |
|---|
| 148 | PRINT *,'Emissivity of bare ground ?' |
|---|
| 149 | call getin("emis",emis) |
|---|
| 150 | write(*,*) " emis = ",emis |
|---|
| 151 | |
|---|
| 152 | pceil=100.0 ! Pascals |
|---|
| 153 | PRINT *,'Ceiling pressure (Pa) ?' |
|---|
| 154 | call getin("pceil",pceil) |
|---|
| 155 | write(*,*) " pceil = ", pceil |
|---|
| 156 | |
|---|
| 157 | mugaz=0. |
|---|
| 158 | cpp=0. |
|---|
| 159 | |
|---|
| 160 | call su_gases |
|---|
| 161 | call calc_cpp_mugaz |
|---|
| 162 | |
|---|
| 163 | |
|---|
| 164 | call inifis(1,llm,0,86400.0,1.0,0.0,0.0,1.0,rad,g,r,cpp) |
|---|
| 165 | |
|---|
| 166 | ! Tracer initialisation |
|---|
| 167 | ! --------------------- |
|---|
| 168 | if (tracer) then |
|---|
| 169 | ! load tracer names from file 'traceur.def' |
|---|
| 170 | open(90,file='traceur.def',status='old',form='formatted',& |
|---|
| 171 | iostat=ierr) |
|---|
| 172 | if (ierr.eq.0) then |
|---|
| 173 | write(*,*) "kcm1d: Reading file traceur.def" |
|---|
| 174 | ! read number of tracers: |
|---|
| 175 | read(90,*,iostat=ierr) nq |
|---|
| 176 | if (ierr.ne.0) then |
|---|
| 177 | write(*,*) "kcm1d: error reading number of tracers" |
|---|
| 178 | write(*,*) " (first line of traceur.def) " |
|---|
| 179 | stop |
|---|
| 180 | else |
|---|
| 181 | ! check that the number of tracers is indeed nqmx |
|---|
| 182 | if (nq.ne.nqmx) then |
|---|
| 183 | write(*,*) "kcm1d: error, wrong number of tracers:" |
|---|
| 184 | write(*,*) "nq=",nq," whereas nqmx=",nqmx |
|---|
| 185 | stop |
|---|
| 186 | endif |
|---|
| 187 | endif |
|---|
| 188 | |
|---|
| 189 | do iq=1,nq |
|---|
| 190 | ! minimal version, just read in the tracer names, 1 per line |
|---|
| 191 | read(90,*,iostat=ierr) tnom(iq) |
|---|
| 192 | if (ierr.ne.0) then |
|---|
| 193 | write(*,*) 'kcm1d: error reading tracer names...' |
|---|
| 194 | stop |
|---|
| 195 | endif |
|---|
| 196 | enddo !of do iq=1,nq |
|---|
| 197 | endif |
|---|
| 198 | !endif |
|---|
| 199 | |
|---|
| 200 | call initracer() |
|---|
| 201 | |
|---|
| 202 | endif |
|---|
| 203 | |
|---|
| 204 | |
|---|
| 205 | do iq=1,nqmx |
|---|
| 206 | do ilay=1,nlayer |
|---|
| 207 | q(ilay,iq) = 0. |
|---|
| 208 | enddo |
|---|
| 209 | enddo |
|---|
| 210 | |
|---|
| 211 | do iq=1,nqmx |
|---|
| 212 | qsurf(iq) = 0. |
|---|
| 213 | enddo |
|---|
| 214 | |
|---|
| 215 | firstcall = .true. |
|---|
| 216 | lastcall = .false. |
|---|
| 217 | |
|---|
| 218 | iter = 1 |
|---|
| 219 | Tstrat = 60.0 |
|---|
| 220 | dTstrat = 1000.0 |
|---|
| 221 | |
|---|
| 222 | ! --------- |
|---|
| 223 | ! Run model |
|---|
| 224 | ! --------- |
|---|
| 225 | do |
|---|
| 226 | psurf = psurf_n |
|---|
| 227 | |
|---|
| 228 | ! Create vertical profiles |
|---|
| 229 | call kcmprof_fn(psurf,qsurf(1),tsurf, & |
|---|
| 230 | tstrat,play,plev,zlay,temp,q(:,1),muvar(1)) |
|---|
| 231 | |
|---|
| 232 | ! Run radiative transfer |
|---|
| 233 | call callcorrk(ngridmx,nlayer,q,nqmx,qsurf, & |
|---|
| 234 | albedo,emis,mu0,plev,play,temp, & |
|---|
| 235 | tsurf,fract,dist_star,aerosol,cpp,muvar, & |
|---|
| 236 | dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, & |
|---|
| 237 | fluxabs_sw,fluxtop_dn,reffrad,tau_col, & |
|---|
| 238 | cloudfrac,totcloudfrac,.false.,firstcall,lastcall) |
|---|
| 239 | |
|---|
| 240 | ! Iterate stratospheric temperature |
|---|
| 241 | print*,'Tstrat = ',Tstrat |
|---|
| 242 | dTstrat = Tstrat |
|---|
| 243 | !Tstrat = Tsurf*(0.2786*(psurf/100000.)**(-1.123))**0.25 |
|---|
| 244 | ! skin temperature (gray approx.) using analytic pure H2 expression |
|---|
| 245 | !Tstrat = (fluxabs_sw/(2*sigma))**0.25 ! skin temperature (gray approx.) |
|---|
| 246 | Tstrat = (fluxtop_lw/(2*sigma))**0.25 ! skin temperature (gray approx.) |
|---|
| 247 | dTstrat = dTstrat-Tstrat |
|---|
| 248 | |
|---|
| 249 | if(abs(dTstrat).lt.1.0)then |
|---|
| 250 | print*,'dTstrat = ',dTstrat |
|---|
| 251 | exit |
|---|
| 252 | endif |
|---|
| 253 | |
|---|
| 254 | iter=iter+1 |
|---|
| 255 | if(iter.eq.100)then |
|---|
| 256 | print*,'Stratosphere failed to converge after' |
|---|
| 257 | print*,'100 iteration, aborting run.' |
|---|
| 258 | call abort |
|---|
| 259 | endif |
|---|
| 260 | |
|---|
| 261 | end do |
|---|
| 262 | |
|---|
| 263 | ! Calculate total atmospheric energy |
|---|
| 264 | Eatmtot=0.0 |
|---|
| 265 | ! call calcenergy_kcm(tsurf,temp,play,plev,qsurf,& |
|---|
| 266 | ! q(:,1),muvar,Eatmtot) |
|---|
| 267 | |
|---|
| 268 | ! ------------------------ |
|---|
| 269 | ! Save data to ascii files |
|---|
| 270 | ! ------------------------ |
|---|
| 271 | |
|---|
| 272 | print*,'Saving profiles...' |
|---|
| 273 | open(115,file='profpres.out',form='formatted') |
|---|
| 274 | open(116,file='proftemp.out',form='formatted') |
|---|
| 275 | open(117,file='profztab.out',form='formatted') |
|---|
| 276 | open(118,file='profqvar.out',form='formatted') |
|---|
| 277 | open(119,file='profvvar.out',form='formatted') |
|---|
| 278 | |
|---|
| 279 | write(115,*) psurf |
|---|
| 280 | write(116,*) tsurf |
|---|
| 281 | write(117,*) 0.0 |
|---|
| 282 | write(118,*) qsurf(1) |
|---|
| 283 | write(119,*) qsurf(1)*(muvar(1)/mH2O) |
|---|
| 284 | do ilay=1,nlayer |
|---|
| 285 | vmr(ilay,1) = q(ilay,1)*(muvar(ilay+1)/mH2O) |
|---|
| 286 | write(115,*) play(ilay) |
|---|
| 287 | write(116,*) temp(ilay) |
|---|
| 288 | write(117,*) zlay(ilay) |
|---|
| 289 | write(118,*) q(ilay,1) |
|---|
| 290 | write(119,*) vmr(ilay,1) |
|---|
| 291 | enddo |
|---|
| 292 | close(115) |
|---|
| 293 | close(116) |
|---|
| 294 | close(117) |
|---|
| 295 | close(118) |
|---|
| 296 | close(119) |
|---|
| 297 | |
|---|
| 298 | print*,'Saving scalars...' |
|---|
| 299 | open(116,file='surf_vals.out') |
|---|
| 300 | write(116,*) tsurf,psurf,fluxtop_dn, & |
|---|
| 301 | fluxabs_sw,fluxtop_lw |
|---|
| 302 | close(116) |
|---|
| 303 | open(111,file='ene_vals.out') |
|---|
| 304 | write(111,*) tsurf,psurf,Eatmtot,Tstrat |
|---|
| 305 | close(111) |
|---|
| 306 | |
|---|
| 307 | print*,'Saving spectra...' |
|---|
| 308 | open(117,file='OLRnu.out') |
|---|
| 309 | do iw=1,L_NSPECTI |
|---|
| 310 | write(117,*) OLR_nu(iw) |
|---|
| 311 | enddo |
|---|
| 312 | close(117) |
|---|
| 313 | |
|---|
| 314 | open(127,file='GSRnu.out') |
|---|
| 315 | do iw=1,L_NSPECTV |
|---|
| 316 | write(127,*) GSR_nu(iw) |
|---|
| 317 | enddo |
|---|
| 318 | close(127) |
|---|
| 319 | |
|---|
| 320 | end program kcm1d |
|---|
