| 1 | PROGRAM radia1d |
|---|
| 2 | USE dimphy |
|---|
| 3 | USE comsaison |
|---|
| 4 | USE comgeomfi |
|---|
| 5 | USE phys_const |
|---|
| 6 | IMPLICIT NONE |
|---|
| 7 | |
|---|
| 8 | #ifdef TOTO |
|---|
| 9 | |
|---|
| 10 | #include "dimensions.h" |
|---|
| 11 | #include "surfdat.h" |
|---|
| 12 | #include "callkeys.h" |
|---|
| 13 | #include "planete.h" |
|---|
| 14 | |
|---|
| 15 | c Arguments : |
|---|
| 16 | c ----------- |
|---|
| 17 | |
|---|
| 18 | INTEGER ngrid,nlayer,nq |
|---|
| 19 | |
|---|
| 20 | REAL plev(nlayermx+1),play(nlayermx) |
|---|
| 21 | REAL temp(nlayermx) |
|---|
| 22 | REAL psrf |
|---|
| 23 | |
|---|
| 24 | c dynamial tendencies |
|---|
| 25 | |
|---|
| 26 | INTEGER l,ierr,aslun,nlevel,iaer |
|---|
| 27 | c |
|---|
| 28 | REAL mu0,fract |
|---|
| 29 | REAL day_ini,time,longitude,latitude |
|---|
| 30 | REAL zlay(nlayermx) |
|---|
| 31 | REAL ztlev(nlayermx+1) |
|---|
| 32 | REAL zplanck |
|---|
| 33 | REAL zrad(nlayermx),zradc(nlayermx) |
|---|
| 34 | REAL zdum1(nlayermx) |
|---|
| 35 | REAL zdum2(nlayermx) |
|---|
| 36 | REAL zdum3(nlayermx) |
|---|
| 37 | REAL zdum4(nlayermx) |
|---|
| 38 | REAL zdum5(nlayermx) |
|---|
| 39 | REAL stephan |
|---|
| 40 | REAL ztim1,ztim2,ztim3 |
|---|
| 41 | REAL zco2,zp |
|---|
| 42 | REAL ls,zmax,tau,tau_tot |
|---|
| 43 | REAL dtlw(nlayermx),dtsw(nlayermx) |
|---|
| 44 | REAL dtlwcl(nlayermx),dtswcl(nlayermx) |
|---|
| 45 | REAL zflux(6) |
|---|
| 46 | |
|---|
| 47 | REAL aerosol(nlayermx,5),cst_aer,pview |
|---|
| 48 | REAL tsurf,tsoil(nsoilmx) |
|---|
| 49 | REAL co2ice,albedo(2),emis |
|---|
| 50 | REAL dtrad(nlayermx),fluxrad |
|---|
| 51 | |
|---|
| 52 | DATA stephan/5.67e-08/ |
|---|
| 53 | DATA psrf,zmax/775.,9./ |
|---|
| 54 | DATA ls,time,latitude,longitude/0.,12.,0.,0./ |
|---|
| 55 | DATA diurnal/.true./ |
|---|
| 56 | DATA tau/.5/ |
|---|
| 57 | |
|---|
| 58 | c WARNING declin and dist_sol are prescribed instead of Ls |
|---|
| 59 | DATA declin,dist_sol/-24.8,1.4/ |
|---|
| 60 | |
|---|
| 61 | c----------------------------------------------------------------------- |
|---|
| 62 | c 1. Initialisations : |
|---|
| 63 | c -------------------- |
|---|
| 64 | |
|---|
| 65 | PRINT*,'tau?' |
|---|
| 66 | READ(5,*) tau |
|---|
| 67 | PRINT*,'time 0 to 24 h' |
|---|
| 68 | READ(5,*) time |
|---|
| 69 | PRINT*,'latitude (degrees)' |
|---|
| 70 | READ(5,*) latitude |
|---|
| 71 | |
|---|
| 72 | pi=2.*asin(1.) |
|---|
| 73 | ls=ls*pi/180. |
|---|
| 74 | time=time/24. |
|---|
| 75 | latitude=latitude*pi/180. |
|---|
| 76 | longitude=longitude*pi/180. |
|---|
| 77 | declin=declin*pi/180. |
|---|
| 78 | |
|---|
| 79 | ngrid=1 |
|---|
| 80 | nlayer=nlayermx |
|---|
| 81 | nlevel=nlayer+1 |
|---|
| 82 | |
|---|
| 83 | rad=3397200. |
|---|
| 84 | omeg=4.*asin(1.)/(88775.) |
|---|
| 85 | g=3.72 |
|---|
| 86 | mugaz=43.49 |
|---|
| 87 | rcp=.256793 |
|---|
| 88 | unjours=88775. |
|---|
| 89 | r = 8.314511*1000./mugaz |
|---|
| 90 | cpp = r/rcp |
|---|
| 91 | PRINT*,'Cp = ',cpp |
|---|
| 92 | PRINT*,'R = ',r |
|---|
| 93 | |
|---|
| 94 | CALL paramdef |
|---|
| 95 | |
|---|
| 96 | tsurf=200. |
|---|
| 97 | DO l=1,nlayer |
|---|
| 98 | zlay(l)=zmax*(l-.5)/nlayer |
|---|
| 99 | plev(l)=psrf*exp(-zmax*(l-1.)/nlayer) |
|---|
| 100 | play(l)=psrf*exp(-zlay(l)) |
|---|
| 101 | ENDDO |
|---|
| 102 | plev(nlevel)=0. |
|---|
| 103 | |
|---|
| 104 | DO l=1,nlayer |
|---|
| 105 | temp(l)=200. |
|---|
| 106 | ENDDO |
|---|
| 107 | |
|---|
| 108 | DO l=1,nlayer |
|---|
| 109 | dtrad(l)=0. |
|---|
| 110 | ENDDO |
|---|
| 111 | fluxrad=0. |
|---|
| 112 | |
|---|
| 113 | albedo(1)=.24 |
|---|
| 114 | albedo(2)=.24 |
|---|
| 115 | emis=1. |
|---|
| 116 | CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) |
|---|
| 117 | |
|---|
| 118 | c a total otempical detemph of .2 for Ps=700Pa |
|---|
| 119 | cst_aer=tau/700. |
|---|
| 120 | pview=1.66 |
|---|
| 121 | CALL sucst(66,19900101,8000000,g,rad,mugaz,rcp,unjours) |
|---|
| 122 | CALL surad(6) |
|---|
| 123 | |
|---|
| 124 | IF (ngrid.NE.ngridmax) THEN |
|---|
| 125 | PRINT*,'STOP in inifis' |
|---|
| 126 | PRINT*,'Probleme de dimenesions :' |
|---|
| 127 | PRINT*,'ngrid = ',ngrid |
|---|
| 128 | PRINT*,'ngridmax = ',ngridmax |
|---|
| 129 | STOP |
|---|
| 130 | ENDIF |
|---|
| 131 | |
|---|
| 132 | c----------------------------------------------------------------------- |
|---|
| 133 | c 2. Calcul of the radiative tendencies : |
|---|
| 134 | c --------------------------------------- |
|---|
| 135 | |
|---|
| 136 | |
|---|
| 137 | c CALL orbite(ls,dist_sol,declin) |
|---|
| 138 | IF(diurnal) THEN |
|---|
| 139 | ztim1=SIN(declin) |
|---|
| 140 | ztim2=COS(declin)*COS(2.*pi*(time-.5)) |
|---|
| 141 | ztim3=-COS(declin)*SIN(2.*pi*(time-.5)) |
|---|
| 142 | CALL solang(ngrid,sin(longitude),cos(longitude), |
|---|
| 143 | s sin(latitude),cos(latitude), |
|---|
| 144 | s ztim1,ztim2,ztim3, |
|---|
| 145 | s mu0,fract) |
|---|
| 146 | PRINT*,'time, declin, sinlon,coslon,sinlat,coslat' |
|---|
| 147 | PRINT*,time, declin,sin(longitude),cos(longitude), |
|---|
| 148 | s sin(latitude),cos(latitude) |
|---|
| 149 | ELSE |
|---|
| 150 | CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad) |
|---|
| 151 | ENDIF |
|---|
| 152 | |
|---|
| 153 | c 2.2 Calcul of the radiative tendencies and fluxes: |
|---|
| 154 | c -------------------------------------------------- |
|---|
| 155 | |
|---|
| 156 | c 2.1.1 layers |
|---|
| 157 | tau_tot=0. |
|---|
| 158 | DO l=1,nlayer |
|---|
| 159 | zp=plev(1)/play(l) |
|---|
| 160 | aerosol(l,1)=amax1(cst_aer* |
|---|
| 161 | s (plev(l)-plev(l+1)) |
|---|
| 162 | s *exp(.03*(1.-zp)),1.e-33) |
|---|
| 163 | aerosol(l,1)=amax1( |
|---|
| 164 | s tau*(plev(l)-plev(l+1))/plev(1) |
|---|
| 165 | s *exp(.007*(1.-zp)),1.e-33) |
|---|
| 166 | tau_tot=tau_tot+aerosol(l,1) |
|---|
| 167 | ENDDO |
|---|
| 168 | PRINT*,'tau total',tau_tot |
|---|
| 169 | DO l=1,nlayermx |
|---|
| 170 | aerosol(l,1)=aerosol(l,1)*tau/tau_tot |
|---|
| 171 | ENDDO |
|---|
| 172 | |
|---|
| 173 | c 2.1.2 levels |
|---|
| 174 | |
|---|
| 175 | c Extrapolation for the air temperature above the surface |
|---|
| 176 | ztlev(1)=temp(1)+ |
|---|
| 177 | s (plev(1)-play(1))* |
|---|
| 178 | s (temp(1)-temp(2))/(play(1)-play(2)) |
|---|
| 179 | |
|---|
| 180 | DO l=2,nlevel-1 |
|---|
| 181 | ztlev(l)=0.5*(temp(l-1)+temp(l)) |
|---|
| 182 | ENDDO |
|---|
| 183 | |
|---|
| 184 | ztlev(nlevel)=temp(nlayer) |
|---|
| 185 | PRINT*,'temp' |
|---|
| 186 | |
|---|
| 187 | DO l=1,nlayer |
|---|
| 188 | zdum1(l)=0. |
|---|
| 189 | zdum2(l)=0. |
|---|
| 190 | zdum3(l)=0. |
|---|
| 191 | zdum4(l)=0. |
|---|
| 192 | zdum5(l)=0. |
|---|
| 193 | DO iaer=2,5 |
|---|
| 194 | aerosol(l,iaer)=0. |
|---|
| 195 | ENDDO |
|---|
| 196 | ENDDO |
|---|
| 197 | |
|---|
| 198 | c 2.3 Calcul of the radiative tendencies and fluxes: |
|---|
| 199 | c -------------------------------------------------- |
|---|
| 200 | |
|---|
| 201 | zco2=.95 |
|---|
| 202 | CALL radite(ngrid,nlayer,0,6,1,1, |
|---|
| 203 | $ aerosol,albedo,zco2,zdum4,emis, |
|---|
| 204 | $ mu0,zdum5, |
|---|
| 205 | $ plev,play,zdum1,zdum2,zdum3,ztlev,tsurf,temp,pview, |
|---|
| 206 | $ dtlw,dtsw,dtlwcl,dtswcl,zflux,zrad,zradc, |
|---|
| 207 | $ fract,dist_sol) |
|---|
| 208 | PRINT*,'radite' |
|---|
| 209 | |
|---|
| 210 | c 2.4 total flux and tendencies: |
|---|
| 211 | c ------------------------------ |
|---|
| 212 | |
|---|
| 213 | c 2.4.1 fluxes |
|---|
| 214 | |
|---|
| 215 | fluxrad=emis*zflux(4) |
|---|
| 216 | $ +zflux(5)*(1.-albedo(1)) |
|---|
| 217 | $ +zflux(6)*(1.-albedo(2)) |
|---|
| 218 | zplanck=tsurf*tsurf |
|---|
| 219 | zplanck=emis* |
|---|
| 220 | $ stephan*zplanck*zplanck |
|---|
| 221 | fluxrad=fluxrad-zplanck |
|---|
| 222 | |
|---|
| 223 | c 2.4.2 temperature tendencies |
|---|
| 224 | |
|---|
| 225 | DO l=1,nlayer |
|---|
| 226 | dtrad(l)=(dtsw(l)+dtlw(l))/unjours |
|---|
| 227 | dtsw(l)=cpp*dtsw(l)/unjours |
|---|
| 228 | dtlw(l)=dtlw(l) |
|---|
| 229 | ENDDO |
|---|
| 230 | |
|---|
| 231 | |
|---|
| 232 | c 2.5 Transformation of the radiative tendencies: |
|---|
| 233 | c ----------------------------------------------- |
|---|
| 234 | |
|---|
| 235 | PRINT*,'Diagnotique for the radaition' |
|---|
| 236 | PRINT*,'albedo, emissiv, mu0,fract,pview,Frad,Planck' |
|---|
| 237 | PRINT*,albedo(1),emis,mu0,fract,pview,fluxrad,zplanck |
|---|
| 238 | PRINT*,'Tlay Tlev Play Plev z aerosol dT/dt SW dT/dt LW (K/day)' |
|---|
| 239 | PRINT*,'unjours',unjours |
|---|
| 240 | DO l=1,nlayer |
|---|
| 241 | WRITE(*,'(2f7.1,2e11.3,f6.1,3e14.5)') |
|---|
| 242 | s temp(l),ztlev(l),play(l),plev(l),zlay(l), |
|---|
| 243 | s g*aerosol(l,1)/(plev(l)-plev(l+1)),dtsw(l),dtlw(l) |
|---|
| 244 | WRITE(55,'(2e15.5)') -dtlw(l),.001*zlay(l)*r*200./g |
|---|
| 245 | WRITE(56,'(2e15.5)') dtsw(l),zlay(l) |
|---|
| 246 | ENDDO |
|---|
| 247 | |
|---|
| 248 | #endif |
|---|
| 249 | |
|---|
| 250 | END |
|---|
