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