Changeset 253 for trunk/LMDZ.GENERIC/libf/phystd
- Timestamp:
- Aug 2, 2011, 11:13:07 AM (13 years ago)
- Location:
- trunk/LMDZ.GENERIC/libf/phystd
- Files:
-
- 23 added
- 12 deleted
- 47 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
r135 r253 1 subroutine aeropacity(ngrid,nlayer,nq,pplev,pq, &2 aerosol,reffrad,QREFvis3d,QREFir3d,tau_col )3 4 use radinc_h, only : naerkind 1 Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, & 2 aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky) 3 4 use radinc_h, only : naerkind, L_TAUMAX 5 5 6 6 implicit none … … 46 46 INTEGER ngrid,nlayer,nq 47 47 48 REAL pplay(ngrid,nlayer) 48 49 REAL pplev(ngrid,nlayer+1) 49 50 REAL pq(ngrid,nlayer,nq) … … 53 54 REAL QREFir3d(ngridmx,nlayermx,naerkind) 54 55 55 REAL tauref(ngrid), tau_col(ngrid) 56 REAL tau_col(ngrid) 57 ! REAL tauref(ngrid), tau_col(ngrid) 58 59 real cloudfrac(ngridmx,nlayermx) 60 real aerosol0 56 61 57 62 INTEGER l,ig,iq,iaer … … 68 73 CHARACTER(LEN=20) :: tracername ! to temporarily store text 69 74 70 71 ! identify tracers 75 ! for fixed dust profiles 76 real topdust, expfactor, zp 77 REAL taudusttmp(ngridmx) ! Temporary dust opacity used before scaling 78 79 ! BENJAMIN MODIFS 80 real CLFtot,CLF1,CLF2 81 real totcloudfrac(ngridmx) 82 logical clearsky 83 84 ! identify tracers 72 85 IF (firstcall) THEN 73 86 74 ! are these tests of any real use ?87 ! are these tests of any real use ? 75 88 IF(ngrid.NE.ngridmx) THEN 76 89 PRINT*,'STOP in aeropacity!' … … 86 99 stop 87 100 endif 88 89 101 90 102 do iq=1,nqmx … … 101 113 write(*,*) " i_h2oice=",i_h2oice 102 114 103 !if((i_h2oice.lt.1) .and. (naerkind.gt.1))then 104 ! print*,'naerkind > 1 but no h2o ice, this will not work.' 105 ! stop 106 !endif 115 if(watercond.and.(.not.aerofixed))then 116 if(naerkind.lt.2)then 117 print*,'Cannot have active H2O clouds with naerkind less than 2!' 118 call abort 119 endif 120 endif 121 122 if(dusttau.gt.0.0)then 123 if(naerkind.lt.3)then 124 print*,'Cannot have active dust with naerkind less than 3!' 125 call abort 126 endif 127 endif 107 128 108 129 firstcall=.false. … … 110 131 111 132 112 aerosol(:,:,:)=0.0113 114 133 DO iaer = 1, naerkind ! Loop on aerosol kind (NOT tracer) 115 134 ! --------------------------------------------------------- 116 aerkind: SELECT CASE (iaer)117 !================================================================== 118 CASE(1) aerkind! CO2 ice (iaer=1)135 aerkind: SELECT CASE (iaer) 136 !================================================================== 137 CASE(1) aerkind ! CO2 ice (iaer=1) 119 138 !================================================================== 120 139 121 140 ! 1. Initialization 122 CALL zerophys(ngrid*nlayer,aerosol(1,1,iaer))141 aerosol(:,:,iaer)=0.0 123 142 124 143 ! 2. Opacity calculation 125 if (aerofixed.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed 126 do ig=1, ngrid 127 do l=1,nlayer 128 aerosol(ig,l,iaer) =1.e-9 129 end do 130 aerosol(ig,14,iaer)=1.e-9 ! single cloud layer option 131 end do 132 else 133 DO ig=1, ngrid 134 DO l=1,nlayer 135 aerosol(ig,l,iaer) = & 136 ( 0.75 * QREFvis3d(ig,l,iaer) / & 137 ( rho_co2 * reffrad(ig,l,iaer) ) ) * & 138 ( pq(ig,l,i_co2ice) + 1.E-9 ) * & 139 ( pplev(ig,l) - pplev(ig,l+1) ) / g 140 ENDDO 141 ENDDO 142 end if 143 144 !================================================================== 145 CASE(2) aerkind ! Water ice crystals (iaer=2) 146 !================================================================== 147 144 if (aerofixed.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed 145 do ig=1, ngrid 146 do l=1,nlayer 147 aerosol(ig,l,iaer)=1.e-9 148 end do 149 !aerosol(ig,12,iaer)=4.0 ! single cloud layer option 150 end do 151 else 152 DO ig=1, ngrid 153 DO l=1,nlayer-1 ! to stop the rad tran bug 154 155 aerosol0 = & 156 ( 0.75 * QREFvis3d(ig,l,iaer) / & 157 ( rho_co2 * reffrad(ig,l,iaer) ) ) * & 158 ( pq(ig,l,i_co2ice) + 1.E-9 ) * & 159 ( pplev(ig,l) - pplev(ig,l+1) ) / g 160 aerosol0 = max(aerosol0,1.e-9) 161 aerosol0 = min(aerosol0,L_TAUMAX) 162 aerosol(ig,l,iaer) = aerosol0 163 ! aerosol(ig,l,iaer) = 0.0 164 165 ! using cloud fraction 166 ! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF)) 167 ! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX) 168 169 170 ENDDO 171 ENDDO 172 end if 173 174 !================================================================== 175 CASE(2) aerkind ! Water ice / liquid (iaer=2) 176 !================================================================== 148 177 149 178 ! 1. Initialization 150 CALL zerophys(ngrid*nlayer,aerosol(1,1,iaer))179 aerosol(:,:,iaer)=0.0 151 180 152 181 ! 2. Opacity calculation 153 ! if (1.eq.1) then 154 if (aerofixed.or.(i_h2oice.eq.0)) then 155 ! print*,'No H2O clouds in the radiative transfer!' 156 do ig=1, ngrid 157 do l=1,nlayer 158 aerosol(ig,l,iaer) =1.e-9 159 end do 160 aerosol(ig,5,iaer)=1.e-9 ! single cloud layer option 161 end do 162 else 163 DO ig=1, ngrid 164 DO l=1,nlayer 165 aerosol(ig,l,iaer) = & 166 ( 0.75 * QREFvis3d(ig,l,iaer) / & 167 ( rho_ice * reffrad(ig,l,iaer) ) ) * & 168 ( pq(ig,l,i_h2oice) + 1.E-9 ) * & 169 ( pplev(ig,l) - pplev(ig,l+1) ) / g 170 ENDDO 171 ENDDO 172 173 end if 174 175 176 !================================================================== 177 END SELECT aerkind 178 182 if (aerofixed.or.(i_h2oice.eq.0).or.clearsky) then 183 do ig=1, ngrid ! to stop the rad tran bug 184 do l=1,nlayer 185 aerosol(ig,l,iaer) =1.e-9 186 end do 187 end do 188 else 189 do ig=1, ngrid 190 do l=1,nlayer-1 ! to stop the rad tran bug 191 192 if(CLFvarying)then 193 CLF1 = max(cloudfrac(ig,l),0.01) 194 CLFtot = max(totcloudfrac(ig),0.01) 195 CLF2 = CLF1/CLFtot 196 CLF2 = min(1.,CLF2) 197 else 198 CLF1=1.0 199 CLF2=CLFfixval 200 endif 201 202 aerosol0 = & 203 ( 0.75 * QREFvis3d(ig,l,iaer) / & 204 ( rho_ice * reffrad(ig,l,iaer) ) ) * & 205 ( pq(ig,l,i_h2oice) + 1.E-9 ) * & 206 ( pplev(ig,l) - pplev(ig,l+1) ) / g / & 207 CLF1 208 209 aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0)) 210 ! why no division in exponential? 211 ! because its already done in aerosol0 212 213 aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9) 214 aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0) 215 216 enddo 217 enddo 218 end if 219 220 221 !================================================================== 222 CASE(3) aerkind ! Dust (iaer=3) 223 !================================================================== 224 225 ! 1. Initialization 226 aerosol(:,:,iaer)=0.0 227 228 topdust=10.0 ! km 229 230 print*,'WARNING, dust is experimental for Early Mars' 231 232 ! 2. Opacity calculation 233 234 do l=1,nlayer 235 do ig=1,ngrid-1 ! to stop the rad tran bug 236 zp=700./pplay(ig,l) 237 aerosol(ig,l,1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) & 238 *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 ) & 239 *QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer) 240 ! takes into account particle variation 241 enddo 242 enddo 243 244 !================================================================== 245 END SELECT aerkind 179 246 ENDDO ! iaer (loop on aerosol kind) 180 247 … … 183 250 ! Column integrated visible optical depth in each point (used for diagnostic) 184 251 185 call zerophys(ngrid,tau_col)252 tau_col(:)=0.0 186 253 do iaer = 1, naerkind 187 254 do l=1,nlayer … … 192 259 end do 193 260 261 do ig=1, ngrid 262 do l=1,nlayer 263 do iaer = 1, naerkind 264 if(aerosol(ig,l,iaer).gt.1.e3)then 265 print*,'WARNING: aerosol=',aerosol(ig,l,iaer) 266 print*,'at ig=',ig,', l=',l,', iaer=',iaer 267 print*,'QREFvis3d=',QREFvis3d(ig,l,iaer) 268 print*,'reffrad=',reffrad(ig,l,iaer) 269 endif 270 end do 271 end do 272 end do 273 274 do ig=1, ngrid 275 if(tau_col(ig).gt.1.e3)then 276 print*,'WARNING: tau_col=',tau_col(ig) 277 print*,'at ig=',ig 278 print*,'aerosol=',aerosol(ig,:,:) 279 print*,'QREFvis3d=',QREFvis3d(ig,:,:) 280 print*,'reffrad=',reffrad(ig,:,:) 281 endif 282 end do 283 194 284 return 195 285 end subroutine aeropacity -
trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90
r135 r253 1 subroutine aeroptproperties(ngrid,nlayer,reffrad,nueffrad, & 2 QVISsQREF3d,omegaVIS3d,gVIS3d, & 3 QIRsQREF3d,omegaIR3d,gIR3d, & 4 QREFvis3d,QREFir3d) 1 SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad, & 2 QVISsQREF3d,omegaVIS3d,gVIS3d, & 3 QIRsQREF3d,omegaIR3d,gIR3d, & 4 QREFvis3d,QREFir3d)!, & 5 ! omegaREFvis3d,omegaREFir3d) 5 6 6 7 use radinc_h, only: L_NSPECTI,L_NSPECTV,naerkind,nsizemax 7 8 use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir 8 use radcommon_h, only: radiustab,nsize,qrefvis,qrefir 9 use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir 10 use radcommon_h, only: radiustab,nsize 9 11 10 12 implicit none 11 13 12 !================================================================== 13 ! 14 ! Purpose 15 ! ------- 16 ! Compute the scattering parameters in each grid 17 ! box, depending on aerosol grain sizes. 14 ! ============================================================= 15 ! Aerosol Optical Properties 18 16 ! 19 ! Notes 20 ! ----- 21 ! Don't forget to set the value of varyingnueff below; If 22 ! the effective variance of the distribution for the given 23 ! aerosol is considered homogeneous in the atmosphere, please 24 ! set varyingnueff(iaer) to .false. Resulting computational 25 ! time will be much better. 26 ! 27 ! Authors 28 ! ------- 29 ! J.-B. Madeleine, F. Montmessin 17 ! Description: 18 ! Compute the scattering parameters in each grid 19 ! box, depending on aerosol grain sizes. Log-normal size 20 ! distribution and Gauss-Legendre integration are used. 21 22 ! Parameters: 23 ! Don't forget to set the value of varyingnueff below; If 24 ! the effective variance of the distribution for the given 25 ! aerosol is considered homogeneous in the atmosphere, please 26 ! set varyingnueff(iaer) to .false. Resulting computational 27 ! time will be much better. 28 29 ! Authors: J.-B. Madeleine, F. Forget, F. Montmessin 30 30 ! Slightly modified and converted to F90 by R. Wordsworth (2009) 31 31 ! Varying nueff section removed by R. Wordsworth for simplicity 32 ! 33 !================================================================== 32 ! ============================================================== 34 33 35 34 #include "dimensions.h" … … 40 39 ! --------------- 41 40 41 42 42 43 ! ============================================================= 43 44 LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false. 44 45 ! ============================================================= 45 46 46 ! Radius axis used for integration 47 DOUBLE PRECISION :: radiusint(nsizemax+1) 48 ! Min. and max radii of the interpolation grid (in METERS) 49 REAL, PARAMETER :: refftabmin = 2e-8 47 ! Min. and max radius of the interpolation grid (in METERS) 48 REAL, PARAMETER :: refftabmin = 2e-8 !2e-8 50 49 ! REAL, PARAMETER :: refftabmax = 35e-6 51 REAL, PARAMETER :: refftabmax = 1e-3 ! CHANGED BY RDW50 REAL, PARAMETER :: refftabmax = 1e-3 52 51 ! Log of the min and max variance of the interpolation grid 53 52 REAL, PARAMETER :: nuefftabmin = -4.6 54 53 REAL, PARAMETER :: nuefftabmax = 0. 55 ! Number of effective radi iof the interpolation grid54 ! Number of effective radius of the interpolation grid 56 55 INTEGER, PARAMETER :: refftabsize = 200 57 56 ! Number of effective variances of the interpolation grid 58 57 ! INTEGER, PARAMETER :: nuefftabsize = 100 59 INTEGER, PARAMETER :: nuefftabsize = 1 ! CHANGED BY RDW58 INTEGER, PARAMETER :: nuefftabsize = 1 60 59 ! Interpolation grid indices (reff,nueff) 61 60 INTEGER :: grid_i,grid_j 62 ! Volume ratio of the look-up table (different in VIS and IR)63 DOUBLE PRECISION :: vrat64 ! r_g and sigma_g for the log-normal distribution65 ! as defined by [hansen_1974]66 REAL :: r_g,sigma_g67 ! Error function used for integration68 DOUBLE PRECISION :: derf69 ! Density function f(x)dx of the log-normal distribution70 REAL :: dfi71 DOUBLE PRECISION :: dfi_tmp(nsizemax+1)72 61 ! Intermediate variable 73 62 REAL :: var_tmp,var3d_tmp(ngridmx,nlayermx) 74 63 ! Bilinear interpolation factors 75 64 REAL :: kx,ky,k1,k2,k3,k4 65 ! Size distribution parameters 66 REAL :: sizedistk1,sizedistk2 67 ! Pi! 68 REAL,SAVE :: pi 69 ! Variables used by the Gauss-Legendre integration: 70 INTEGER radius_id,gausind 71 REAL kint 72 REAL drad 73 INTEGER, PARAMETER :: ngau = 10 74 REAL weightgaus(ngau),radgaus(ngau) 75 SAVE weightgaus,radgaus 76 ! DATA weightgaus/.2955242247,.2692667193,.2190863625,.1494513491,.0666713443/ 77 ! DATA radgaus/.1488743389,.4333953941,.6794095682,.8650633666,.9739065285/ 78 DATA radgaus/0.07652652113350,0.22778585114165, & 79 0.37370608871528,0.51086700195146, & 80 0.63605368072468,0.74633190646476, & 81 0.83911697181213,0.91223442826796, & 82 0.96397192726078,0.99312859919241/ 83 84 DATA weightgaus/0.15275338723120,0.14917298659407, & 85 0.14209610937519,0.13168863843930, & 86 0.11819453196154,0.10193011980823, & 87 0.08327674160932,0.06267204829828, & 88 0.04060142982019,0.01761400714091/ 76 89 ! Indices 77 90 INTEGER :: i,j,k,l,m,iaer,idomain … … 82 95 83 96 ! Radius axis of the interpolation grid 84 DOUBLE PRECISION,SAVE :: refftab(refftabsize)97 REAL,SAVE :: refftab(refftabsize) 85 98 ! Variance axis of the interpolation grid 86 DOUBLE PRECISION,SAVE :: nuefftab(nuefftabsize)99 REAL,SAVE :: nuefftab(nuefftabsize) 87 100 ! Volume ratio of the grid 88 DOUBLE PRECISION,SAVE :: logvratgrid,vratgrid101 REAL,SAVE :: logvratgrid,vratgrid 89 102 ! Grid used to remember which calculation is done 90 LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2) & 91 = .false. 103 LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2) = .false. 92 104 ! Optical properties of the grid (VISIBLE) 93 REAL,SAVE :: epVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) 105 REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) 106 REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) 107 REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) 94 108 REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) 95 109 REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind) 96 110 ! Optical properties of the grid (INFRARED) 97 REAL,SAVE :: epIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) 111 REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) 112 REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) 113 REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) 98 114 REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) 99 115 REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind) 100 116 ! Optical properties of the grid (REFERENCE WAVELENGTHS) 101 REAL,SAVE :: eprefVISgrid(refftabsize,nuefftabsize,naerkind) 102 REAL,SAVE :: eprefIRgrid(refftabsize,nuefftabsize,naerkind) 117 REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind) 118 REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind) 119 REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind) 120 REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind) 121 REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind) 122 REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind) 103 123 ! Firstcall 104 124 LOGICAL,SAVE :: firstcall = .true. 125 ! Variables used by the Gauss-Legendre integration: 126 REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2) 127 REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau) 128 REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau) 129 130 REAL,SAVE :: radGAUSa(ngau,naerkind,2) 131 REAL,SAVE :: radGAUSb(ngau,naerkind,2) 132 133 REAL,SAVE :: qsqrefVISa(L_NSPECTV,ngau,naerkind) 134 REAL,SAVE :: qrefVISa(ngau,naerkind) 135 REAL,SAVE :: qsqrefVISb(L_NSPECTV,ngau,naerkind) 136 REAL,SAVE :: qrefVISb(ngau,naerkind) 137 REAL,SAVE :: omegVISa(L_NSPECTV,ngau,naerkind) 138 REAL,SAVE :: omegrefVISa(ngau,naerkind) 139 REAL,SAVE :: omegVISb(L_NSPECTV,ngau,naerkind) 140 REAL,SAVE :: omegrefVISb(ngau,naerkind) 141 REAL,SAVE :: gVISa(L_NSPECTV,ngau,naerkind) 142 REAL,SAVE :: gVISb(L_NSPECTV,ngau,naerkind) 143 144 REAL,SAVE :: qsqrefIRa(L_NSPECTI,ngau,naerkind) 145 REAL,SAVE :: qrefIRa(ngau,naerkind) 146 REAL,SAVE :: qsqrefIRb(L_NSPECTI,ngau,naerkind) 147 REAL,SAVE :: qrefIRb(ngau,naerkind) 148 REAL,SAVE :: omegIRa(L_NSPECTI,ngau,naerkind) 149 REAL,SAVE :: omegrefIRa(ngau,naerkind) 150 REAL,SAVE :: omegIRb(L_NSPECTI,ngau,naerkind) 151 REAL,SAVE :: omegrefIRb(ngau,naerkind) 152 REAL,SAVE :: gIRa(L_NSPECTI,ngau,naerkind) 153 REAL,SAVE :: gIRb(L_NSPECTI,ngau,naerkind) 154 155 REAL :: radiusm 156 REAL :: radiusr 105 157 106 158 ! Inputs … … 127 179 REAL :: QREFir3d(ngridmx,nlayermx,naerkind) 128 180 181 REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind) 182 REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) 183 129 184 DO iaer = 1, naerkind ! Loop on aerosol kind 130 131 185 IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN 132 133 186 !================================================================== 134 187 ! If there is one single particle size, optical … … 148 201 QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1) 149 202 QREFir3d(ig,lg,iaer)=QREFir(iaer,1) 203 omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1) 204 omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1) 150 205 ENDDO 151 206 ENDDO 207 152 208 153 209 if (firstcall) then … … 162 218 ! -------------------------------------------------- 163 219 IF (firstcall) THEN 164 ! 1.1 Effective radius 220 221 ! 1.1 Pi! 222 pi = 2. * asin(1.e0) 223 224 ! 1.2 Effective radius 165 225 refftab(1) = refftabmin 166 226 refftab(refftabsize) = refftabmax 167 227 168 logvratgrid = log(refftabmax/refftabmin) / & 169 float(refftabsize-1)*3. 228 logvratgrid = log(refftabmax/refftabmin) / float(refftabsize-1)*3. 170 229 vratgrid = exp(logvratgrid) 171 230 … … 174 233 enddo 175 234 176 ! 1.2 Effective variance 177 do i = 0, nuefftabsize-1 178 nuefftab(i+1) = exp( nuefftabmin + & 179 i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) ) 180 enddo 235 ! 1.3 Effective variance 236 if(nuefftabsize.eq.1)then ! addded by RDW 237 print*,'Warning: no variance range in aeroptproperties' 238 nuefftab(1)=0.2 239 else 240 do i = 0, nuefftabsize-1 241 nuefftab(i+1) = exp( nuefftabmin + i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) ) 242 enddo 243 endif 244 181 245 firstcall = .false. 182 246 ENDIF 183 247 184 248 ! 1.4 Radius middle point and range for Gauss integration 249 radiusm=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + radiustab(iaer,idomain,1)) 250 radiusr=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - radiustab(iaer,idomain,1)) 251 252 ! 1.5 Interpolating data at the Gauss quadrature points: 253 DO gausind=1,ngau 254 drad=radiusr*radgaus(gausind) 255 radGAUSa(gausind,iaer,idomain)=radiusm-drad 256 257 radius_id=minloc(abs(radiustab(iaer,idomain,:) - (radiusm-drad)),DIM=1) 258 IF ((radiustab(iaer,idomain,radius_id) - (radiusm-drad)).GT.0) THEN 259 radius_id=radius_id-1 260 ENDIF 261 IF (radius_id.GE.nsize(iaer,idomain)) THEN 262 radius_id=nsize(iaer,idomain)-1 263 kint = 1. 264 ELSEIF (radius_id.LT.1) THEN 265 radius_id=1 266 kint = 0. 267 ELSE 268 kint = ( (radiusm-drad) - & 269 radiustab(iaer,idomain,radius_id) ) / & 270 ( radiustab(iaer,idomain,radius_id+1) - & 271 radiustab(iaer,idomain,radius_id) ) 272 ENDIF 273 IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------------- 274 DO m=1,L_NSPECTV 275 qsqrefVISa(m,gausind,iaer)= & 276 (1-kint)*QVISsQREF(m,iaer,radius_id) + & 277 kint*QVISsQREF(m,iaer,radius_id+1) 278 omegVISa(m,gausind,iaer)= & 279 (1-kint)*omegaVIS(m,iaer,radius_id) + & 280 kint*omegaVIS(m,iaer,radius_id+1) 281 gVISa(m,gausind,iaer)= & 282 (1-kint)*gVIS(m,iaer,radius_id) + & 283 kint*gVIS(m,iaer,radius_id+1) 284 ENDDO 285 qrefVISa(gausind,iaer)= & 286 (1-kint)*QREFvis(iaer,radius_id) + & 287 kint*QREFvis(iaer,radius_id+1) 288 omegrefVISa(gausind,iaer)= & 289 (1-kint)*omegaREFvis(iaer,radius_id) + & 290 kint*omegaREFvis(iaer,radius_id+1) 291 ELSE ! INFRARED DOMAIN ---------------------------------- 292 DO m=1,L_NSPECTI 293 qsqrefIRa(m,gausind,iaer)= & 294 (1-kint)*QIRsQREF(m,iaer,radius_id) + & 295 kint*QIRsQREF(m,iaer,radius_id+1) 296 omegIRa(m,gausind,iaer)= & 297 (1-kint)*omegaIR(m,iaer,radius_id) + & 298 kint*omegaIR(m,iaer,radius_id+1) 299 gIRa(m,gausind,iaer)= & 300 (1-kint)*gIR(m,iaer,radius_id) + & 301 kint*gIR(m,iaer,radius_id+1) 302 ENDDO 303 qrefIRa(gausind,iaer)= & 304 (1-kint)*QREFir(iaer,radius_id) + & 305 kint*QREFir(iaer,radius_id+1) 306 omegrefIRa(gausind,iaer)= & 307 (1-kint)*omegaREFir(iaer,radius_id) + & 308 kint*omegaREFir(iaer,radius_id+1) 309 ENDIF 310 ENDDO 311 312 DO gausind=1,ngau 313 drad=radiusr*radgaus(gausind) 314 radGAUSb(gausind,iaer,idomain)=radiusm+drad 315 316 radius_id=minloc(abs(radiustab(iaer,idomain,:) - & 317 (radiusm+drad)),DIM=1) 318 IF ((radiustab(iaer,idomain,radius_id) - & 319 (radiusm+drad)).GT.0) THEN 320 radius_id=radius_id-1 321 ENDIF 322 IF (radius_id.GE.nsize(iaer,idomain)) THEN 323 radius_id=nsize(iaer,idomain)-1 324 kint = 1. 325 ELSEIF (radius_id.LT.1) THEN 326 radius_id=1 327 kint = 0. 328 ELSE 329 kint = ( (radiusm+drad) - & 330 radiustab(iaer,idomain,radius_id) ) / & 331 ( radiustab(iaer,idomain,radius_id+1) - & 332 radiustab(iaer,idomain,radius_id) ) 333 ENDIF 334 IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------------- 335 DO m=1,L_NSPECTV 336 qsqrefVISb(m,gausind,iaer)= & 337 (1-kint)*QVISsQREF(m,iaer,radius_id) + & 338 kint*QVISsQREF(m,iaer,radius_id+1) 339 omegVISb(m,gausind,iaer)= & 340 (1-kint)*omegaVIS(m,iaer,radius_id) + & 341 kint*omegaVIS(m,iaer,radius_id+1) 342 gVISb(m,gausind,iaer)= & 343 (1-kint)*gVIS(m,iaer,radius_id) + & 344 kint*gVIS(m,iaer,radius_id+1) 345 ENDDO 346 qrefVISb(gausind,iaer)= & 347 (1-kint)*QREFvis(iaer,radius_id) + & 348 kint*QREFvis(iaer,radius_id+1) 349 omegrefVISb(gausind,iaer)= & 350 (1-kint)*omegaREFvis(iaer,radius_id) + & 351 kint*omegaREFvis(iaer,radius_id+1) 352 ELSE ! INFRARED DOMAIN ---------------------------------- 353 DO m=1,L_NSPECTI 354 qsqrefIRb(m,gausind,iaer)= & 355 (1-kint)*QIRsQREF(m,iaer,radius_id) + & 356 kint*QIRsQREF(m,iaer,radius_id+1) 357 omegIRb(m,gausind,iaer)= & 358 (1-kint)*omegaIR(m,iaer,radius_id) + & 359 kint*omegaIR(m,iaer,radius_id+1) 360 gIRb(m,gausind,iaer)= & 361 (1-kint)*gIR(m,iaer,radius_id) + & 362 kint*gIR(m,iaer,radius_id+1) 363 ENDDO 364 qrefIRb(gausind,iaer)= & 365 (1-kint)*QREFir(iaer,radius_id) + & 366 kint*QREFir(iaer,radius_id+1) 367 omegrefIRb(gausind,iaer)= & 368 (1-kint)*omegaREFir(iaer,radius_id) + & 369 kint*omegaREFir(iaer,radius_id+1) 370 ENDIF 371 ENDDO 372 373 !================================================================== 374 ! CONSTANT NUEFF FROM HERE 375 !================================================================== 185 376 186 377 ! 2. Compute the scattering parameters using linear 187 378 ! interpolation over grain sizes and constant nueff 188 379 ! --------------------------------------------------- 189 190 ! 2.1 Initialization191 192 vrat = log(radiustab(iaer,idomain,nsize(iaer,idomain)) / &193 radiustab(iaer,idomain,1)) / &194 float(nsize(iaer,idomain)-1)*3.195 vrat = exp(vrat)196 197 radiusint(1) = 1.e-9198 DO i = 2,nsize(iaer,idomain)199 radiusint(i) = ( (2.*vrat) / (vrat+1.) )**(1./3.) * &200 radiustab(iaer,idomain,i-1)201 ENDDO202 radiusint(nsize(iaer,idomain)+1) = 1.e-2203 380 204 381 DO lg = 1,nlayer … … 210 387 grid_i=floor(var_tmp) 211 388 IF (grid_i.GE.refftabsize) THEN 212 WRITE(*,*) 'Warning: Aerosol particle size in grid box #' 213 WRITE(*,*) ig,' is too large to be used by the ' 214 WRITE(*,*) 'radiative transfer; please extend the ' 215 WRITE(*,*) 'interpolation grid to larger sizes.' 216 389 ! WRITE(*,*) 'Warning: particle size in grid box #' 390 ! WRITE(*,*) ig,' is too large to be used by the ' 391 ! WRITE(*,*) 'radiative transfer; please extend the ' 392 ! WRITE(*,*) 'interpolation grid to larger grain sizes.' 217 393 grid_i=refftabsize-1 218 394 kx = 1. 219 395 ELSEIF (grid_i.LT.1) THEN 220 WRITE(*,*) 'Warning: Aerosolparticle size in grid box #'221 WRITE(*,*) ig,' is too small to be used by the '222 WRITE(*,*) 'radiative transfer; please extend the '223 WRITE(*,*) 'interpolation grid to smallersizes.'396 ! WRITE(*,*) 'Warning: particle size in grid box #' 397 ! WRITE(*,*) ig,' is too small to be used by the ' 398 ! WRITE(*,*) 'radiative transfer; please extend the ' 399 ! WRITE(*,*) 'interpolation grid to smaller grain sizes.' 224 400 grid_i=1 225 401 kx = 0. 226 402 ELSE 227 kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) / &228 ( refftab(grid_i+1)-refftab(grid_i) )403 kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) / & 404 ( refftab(grid_i+1)-refftab(grid_i) ) 229 405 ENDIF 230 406 ! 2.3 Integration 231 407 DO j=grid_i,grid_i+1 232 ! 2.3.1 Check if the calculation has been completed408 ! 2.3.1 Check if the calculation has been done 233 409 IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN 234 ! 2.3.2 Compute r_g and sigma_g for the log-normal 235 ! distribution as defined by [hansen_1974], "Light 236 ! scattering in planetary atmospheres", Space 237 ! Science Reviews 16 527-610, p558 238 sigma_g = log(1.+nueffrad(1,1,iaer)) 239 r_g = exp(2.5*sigma_g) 240 sigma_g = sqrt(sigma_g) 241 r_g = refftab(j) / r_g 410 ! 2.3.2 Log-normal dist., r_g and sigma_g are defined 411 ! in [hansen_1974], "Light scattering in planetary 412 ! atmospheres", Space Science Reviews 16 527-610. 413 ! Here, sizedistk1=r_g and sizedistk2=sigma_g^2 414 sizedistk2 = log(1.+nueffrad(1,1,iaer)) 415 sizedistk1 = exp(2.5*sizedistk2) 416 sizedistk1 = refftab(j) / sizedistk1 417 418 normd(j,1,iaer,idomain) = 1e-30 419 DO gausind=1,ngau 420 drad=radiusr*radgaus(gausind) 421 dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1) 422 dista(j,1,iaer,idomain,gausind) = & 423 EXP(-dista(j,1,iaer,idomain,gausind) * & 424 dista(j,1,iaer,idomain,gausind) * & 425 0.5e0/sizedistk2)/(radiusm-drad) 426 dista(j,1,iaer,idomain,gausind) = & 427 dista(j,1,iaer,idomain,gausind) / & 428 (sqrt(2e0*pi*sizedistk2)) 429 430 distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1) 431 distb(j,1,iaer,idomain,gausind) = & 432 EXP(-distb(j,1,iaer,idomain,gausind) * & 433 distb(j,1,iaer,idomain,gausind) * & 434 0.5e0/sizedistk2)/(radiusm+drad) 435 distb(j,1,iaer,idomain,gausind) = & 436 distb(j,1,iaer,idomain,gausind) / & 437 (sqrt(2e0*pi*sizedistk2)) 438 439 normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) + & 440 weightgaus(gausind) * & 441 ( & 442 distb(j,1,iaer,idomain,gausind) * pi * & 443 radGAUSb(gausind,iaer,idomain) * & 444 radGAUSb(gausind,iaer,idomain) + & 445 dista(j,1,iaer,idomain,gausind) * pi * & 446 radGAUSa(gausind,iaer,idomain) * & 447 radGAUSa(gausind,iaer,idomain) & 448 ) 449 ENDDO 242 450 IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------- 243 451 ! 2.3.3.vis Initialization 244 epVISgrid(j,1,:,iaer)=0. 452 qsqrefVISgrid(j,1,:,iaer)=0. 453 qextVISgrid(j,1,:,iaer)=0. 454 qscatVISgrid(j,1,:,iaer)=0. 245 455 omegVISgrid(j,1,:,iaer)=0. 246 456 gVISgrid(j,1,:,iaer)=0. 247 eprefVISgrid(j,1,iaer)=0. 248 ! 2.3.4.vis Log-normal distribution 249 DO l=1,nsize(iaer,idomain)+1 250 dfi_tmp(l) = log(radiusint(l)/r_g) / & 251 sqrt(2.)/sigma_g 457 qrefVISgrid(j,1,iaer)=0. 458 qscatrefVISgrid(j,1,iaer)=0. 459 omegrefVISgrid(j,1,iaer)=0. 460 461 DO gausind=1,ngau 462 DO m=1,L_NSPECTV 463 ! Convolution: 464 qextVISgrid(j,1,m,iaer) = & 465 qextVISgrid(j,1,m,iaer) + & 466 weightgaus(gausind) * & 467 ( & 468 qsqrefVISb(m,gausind,iaer) * & 469 qrefVISb(gausind,iaer) * & 470 pi*radGAUSb(gausind,iaer,idomain) * & 471 radGAUSb(gausind,iaer,idomain) * & 472 distb(j,1,iaer,idomain,gausind) + & 473 qsqrefVISa(m,gausind,iaer) * & 474 qrefVISa(gausind,iaer) * & 475 pi*radGAUSa(gausind,iaer,idomain) * & 476 radGAUSa(gausind,iaer,idomain) * & 477 dista(j,1,iaer,idomain,gausind) & 478 ) 479 qscatVISgrid(j,1,m,iaer) = & 480 qscatVISgrid(j,1,m,iaer) + & 481 weightgaus(gausind) * & 482 ( & 483 omegVISb(m,gausind,iaer) * & 484 qsqrefVISb(m,gausind,iaer) * & 485 qrefVISb(gausind,iaer) * & 486 pi*radGAUSb(gausind,iaer,idomain) * & 487 radGAUSb(gausind,iaer,idomain) * & 488 distb(j,1,iaer,idomain,gausind) + & 489 omegVISa(m,gausind,iaer) * & 490 qsqrefVISa(m,gausind,iaer) * & 491 qrefVISa(gausind,iaer) * & 492 pi*radGAUSa(gausind,iaer,idomain) * & 493 radGAUSa(gausind,iaer,idomain) * & 494 dista(j,1,iaer,idomain,gausind) & 495 ) 496 gVISgrid(j,1,m,iaer) = & 497 gVISgrid(j,1,m,iaer) + & 498 weightgaus(gausind) * & 499 ( & 500 omegVISb(m,gausind,iaer) * & 501 qsqrefVISb(m,gausind,iaer) * & 502 qrefVISb(gausind,iaer) * & 503 gVISb(m,gausind,iaer) * & 504 pi*radGAUSb(gausind,iaer,idomain) * & 505 radGAUSb(gausind,iaer,idomain) * & 506 distb(j,1,iaer,idomain,gausind) + & 507 omegVISa(m,gausind,iaer) * & 508 qsqrefVISa(m,gausind,iaer) * & 509 qrefVISa(gausind,iaer) * & 510 gVISa(m,gausind,iaer) * & 511 pi*radGAUSa(gausind,iaer,idomain) * & 512 radGAUSa(gausind,iaer,idomain) * & 513 dista(j,1,iaer,idomain,gausind) & 514 ) 515 ENDDO 516 qrefVISgrid(j,1,iaer) = & 517 qrefVISgrid(j,1,iaer) + & 518 weightgaus(gausind) * & 519 ( & 520 qrefVISb(gausind,iaer) * & 521 pi*radGAUSb(gausind,iaer,idomain) * & 522 radGAUSb(gausind,iaer,idomain) * & 523 distb(j,1,iaer,idomain,gausind) + & 524 qrefVISa(gausind,iaer) * & 525 pi*radGAUSa(gausind,iaer,idomain) * & 526 radGAUSa(gausind,iaer,idomain) * & 527 dista(j,1,iaer,idomain,gausind) & 528 ) 529 qscatrefVISgrid(j,1,iaer) = & 530 qscatrefVISgrid(j,1,iaer) + & 531 weightgaus(gausind) * & 532 ( & 533 omegrefVISb(gausind,iaer) * & 534 qrefVISb(gausind,iaer) * & 535 pi*radGAUSb(gausind,iaer,idomain) * & 536 radGAUSb(gausind,iaer,idomain) * & 537 distb(j,1,iaer,idomain,gausind) + & 538 omegrefVISa(gausind,iaer) * & 539 qrefVISa(gausind,iaer) * & 540 pi*radGAUSa(gausind,iaer,idomain) * & 541 radGAUSa(gausind,iaer,idomain) * & 542 dista(j,1,iaer,idomain,gausind) & 543 ) 252 544 ENDDO 253 DO l=1,nsize(iaer,idomain) 254 dfi = 0.5*( derf(dfi_tmp(l+1)) - & 255 derf(dfi_tmp(l)) ) 256 DO m=1,L_NSPECTV 257 epVISgrid(j,1,m,iaer) = & 258 epVISgrid(j,1,m,iaer) & 259 + QVISsQREF(m,iaer,l)*dfi 260 omegVISgrid(j,1,m,iaer) = & 261 omegVISgrid(j,1,m,iaer) & 262 + omegaVIS(m,iaer,l)*dfi 263 gVISgrid(j,1,m,iaer) = & 264 gVISgrid(j,1,m,iaer) & 265 + gVIS(m,iaer,l)*dfi 266 ENDDO !L_NSPECTV 267 eprefVISgrid(j,1,iaer) = & 268 eprefVISgrid(j,1,iaer) & 269 + QREFvis(iaer,l)*dfi 270 ENDDO !nsize 545 546 qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) / & 547 normd(j,1,iaer,idomain) 548 qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / & 549 normd(j,1,iaer,idomain) 550 omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / & 551 qrefVISgrid(j,1,iaer) 552 DO m=1,L_NSPECTV 553 qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / & 554 normd(j,1,iaer,idomain) 555 qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / & 556 normd(j,1,iaer,idomain) 557 gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) / & 558 qscatVISgrid(j,1,m,iaer) / & 559 normd(j,1,iaer,idomain) 560 561 qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / & 562 qrefVISgrid(j,1,iaer) 563 omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / & 564 qextVISgrid(j,1,m,iaer) 565 ENDDO 271 566 ELSE ! INFRARED DOMAIN ---------- 272 567 ! 2.3.3.ir Initialization 273 epIRgrid(j,1,:,iaer)=0. 568 qsqrefIRgrid(j,1,:,iaer)=0. 569 qextIRgrid(j,1,:,iaer)=0. 570 qscatIRgrid(j,1,:,iaer)=0. 274 571 omegIRgrid(j,1,:,iaer)=0. 275 572 gIRgrid(j,1,:,iaer)=0. 276 eprefIRgrid(j,1,iaer)=0. 277 ! 2.3.4.ir Log-normal distribution 278 DO l=1,nsize(iaer,idomain)+1 279 dfi_tmp(l) = log(radiusint(l)/r_g) / & 280 sqrt(2.)/sigma_g 573 qrefIRgrid(j,1,iaer)=0. 574 qscatrefIRgrid(j,1,iaer)=0. 575 omegrefIRgrid(j,1,iaer)=0. 576 577 DO gausind=1,ngau 578 DO m=1,L_NSPECTI 579 ! Convolution: 580 qextIRgrid(j,1,m,iaer) = & 581 qextIRgrid(j,1,m,iaer) + & 582 weightgaus(gausind) * & 583 ( & 584 qsqrefIRb(m,gausind,iaer) * & 585 qrefVISb(gausind,iaer) * & 586 pi*radGAUSb(gausind,iaer,idomain) * & 587 radGAUSb(gausind,iaer,idomain) * & 588 distb(j,1,iaer,idomain,gausind) + & 589 qsqrefIRa(m,gausind,iaer) * & 590 qrefVISa(gausind,iaer) * & 591 pi*radGAUSa(gausind,iaer,idomain) * & 592 radGAUSa(gausind,iaer,idomain) * & 593 dista(j,1,iaer,idomain,gausind) & 594 ) 595 qscatIRgrid(j,1,m,iaer) = & 596 qscatIRgrid(j,1,m,iaer) + & 597 weightgaus(gausind) * & 598 ( & 599 omegIRb(m,gausind,iaer) * & 600 qsqrefIRb(m,gausind,iaer) * & 601 qrefVISb(gausind,iaer) * & 602 pi*radGAUSb(gausind,iaer,idomain) * & 603 radGAUSb(gausind,iaer,idomain) * & 604 distb(j,1,iaer,idomain,gausind) + & 605 omegIRa(m,gausind,iaer) * & 606 qsqrefIRa(m,gausind,iaer) * & 607 qrefVISa(gausind,iaer) * & 608 pi*radGAUSa(gausind,iaer,idomain) * & 609 radGAUSa(gausind,iaer,idomain) * & 610 dista(j,1,iaer,idomain,gausind) & 611 ) 612 gIRgrid(j,1,m,iaer) = & 613 gIRgrid(j,1,m,iaer) + & 614 weightgaus(gausind) * & 615 ( & 616 omegIRb(m,gausind,iaer) * & 617 qsqrefIRb(m,gausind,iaer) * & 618 qrefVISb(gausind,iaer) * & 619 gIRb(m,gausind,iaer) * & 620 pi*radGAUSb(gausind,iaer,idomain) * & 621 radGAUSb(gausind,iaer,idomain) * & 622 distb(j,1,iaer,idomain,gausind) + & 623 omegIRa(m,gausind,iaer) * & 624 qsqrefIRa(m,gausind,iaer) * & 625 qrefVISa(gausind,iaer) * & 626 gIRa(m,gausind,iaer) * & 627 pi*radGAUSa(gausind,iaer,idomain) * & 628 radGAUSa(gausind,iaer,idomain) * & 629 dista(j,1,iaer,idomain,gausind) & 630 ) 631 ENDDO 632 qrefIRgrid(j,1,iaer) = & 633 qrefIRgrid(j,1,iaer) + & 634 weightgaus(gausind) * & 635 ( & 636 qrefIRb(gausind,iaer) * & 637 pi*radGAUSb(gausind,iaer,idomain) * & 638 radGAUSb(gausind,iaer,idomain) * & 639 distb(j,1,iaer,idomain,gausind) + & 640 qrefIRa(gausind,iaer) * & 641 pi*radGAUSa(gausind,iaer,idomain) * & 642 radGAUSa(gausind,iaer,idomain) * & 643 dista(j,1,iaer,idomain,gausind) & 644 ) 645 qscatrefIRgrid(j,1,iaer) = & 646 qscatrefIRgrid(j,1,iaer) + & 647 weightgaus(gausind) * & 648 ( & 649 omegrefIRb(gausind,iaer) * & 650 qrefIRb(gausind,iaer) * & 651 pi*radGAUSb(gausind,iaer,idomain) * & 652 radGAUSb(gausind,iaer,idomain) * & 653 distb(j,1,iaer,idomain,gausind) + & 654 omegrefIRa(gausind,iaer) * & 655 qrefIRa(gausind,iaer) * & 656 pi*radGAUSa(gausind,iaer,idomain) * & 657 radGAUSa(gausind,iaer,idomain) * & 658 dista(j,1,iaer,idomain,gausind) & 659 ) 281 660 ENDDO 282 DO l=1,nsize(iaer,idomain) 283 dfi = 0.5*( derf(dfi_tmp(l+1)) - & 284 derf(dfi_tmp(l)) ) 285 DO m=1,L_NSPECTI 286 epIRgrid(j,1,m,iaer) = & 287 epIRgrid(j,1,m,iaer) & 288 + QIRsQREF(m,iaer,l)*dfi 289 omegIRgrid(j,1,m,iaer) = & 290 omegIRgrid(j,1,m,iaer) & 291 + omegaIR(m,iaer,l)*dfi 292 gIRgrid(j,1,m,iaer) = & 293 gIRgrid(j,1,m,iaer) & 294 + gIR(m,iaer,l)*dfi 295 ENDDO !L_NSPECTI 296 eprefIRgrid(j,1,iaer) = & 297 eprefIRgrid(j,1,iaer) & 298 + QREFir(iaer,l)*dfi 299 ENDDO !nsize 661 662 qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) / & 663 normd(j,1,iaer,idomain) 664 qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / & 665 normd(j,1,iaer,idomain) 666 omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / & 667 qrefIRgrid(j,1,iaer) 668 DO m=1,L_NSPECTI 669 qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / & 670 normd(j,1,iaer,idomain) 671 qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / & 672 normd(j,1,iaer,idomain) 673 gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) / & 674 qscatIRgrid(j,1,m,iaer) / & 675 normd(j,1,iaer,idomain) 676 677 qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / & 678 qrefVISgrid(j,1,iaer) 679 omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / & 680 qextIRgrid(j,1,m,iaer) 681 ENDDO 300 682 ENDIF ! -------------------------- 301 683 checkgrid(j,1,iaer,idomain) = .true. … … 307 689 IF (idomain.EQ.1) THEN ! VISIBLE ------------------------ 308 690 DO m=1,L_NSPECTV 309 QVISsQREF3d(ig,lg,m,iaer) =&310 k1*epVISgrid(grid_i,1,m,iaer) +&311 k2*epVISgrid(grid_i+1,1,m,iaer)312 omegaVIS3d(ig,lg,m,iaer) = &313 k1*omegVISgrid(grid_i,1,m,iaer) +&314 k2*omegVISgrid(grid_i+1,1,m,iaer)315 gVIS3d(ig,lg,m,iaer) = &316 k1*gVISgrid(grid_i,1,m,iaer) +&317 k2*gVISgrid(grid_i+1,1,m,iaer)691 QVISsQREF3d(ig,lg,m,iaer) = & 692 k1*qsqrefVISgrid(grid_i,1,m,iaer) + & 693 k2*qsqrefVISgrid(grid_i+1,1,m,iaer) 694 omegaVIS3d(ig,lg,m,iaer) = & 695 k1*omegVISgrid(grid_i,1,m,iaer) + & 696 k2*omegVISgrid(grid_i+1,1,m,iaer) 697 gVIS3d(ig,lg,m,iaer) = & 698 k1*gVISgrid(grid_i,1,m,iaer) + & 699 k2*gVISgrid(grid_i+1,1,m,iaer) 318 700 ENDDO !L_NSPECTV 319 QREFvis3d(ig,lg,iaer) = & 320 k1*eprefVISgrid(grid_i,1,iaer) + & 321 k2*eprefVISgrid(grid_i+1,1,iaer) 701 QREFvis3d(ig,lg,iaer) = & 702 k1*qrefVISgrid(grid_i,1,iaer) + & 703 k2*qrefVISgrid(grid_i+1,1,iaer) 704 omegaREFvis3d(ig,lg,iaer) = & 705 k1*omegrefVISgrid(grid_i,1,iaer) + & 706 k2*omegrefVISgrid(grid_i+1,1,iaer) 322 707 ELSE ! INFRARED ----------------------- 323 708 DO m=1,L_NSPECTI 324 QIRsQREF3d(ig,lg,m,iaer) = &325 k1*epIRgrid(grid_i,1,m,iaer) +&326 k2*epIRgrid(grid_i+1,1,m,iaer)327 omegaIR3d(ig,lg,m,iaer) = &328 k1*omegIRgrid(grid_i,1,m,iaer) +&329 k2*omegIRgrid(grid_i+1,1,m,iaer)330 gIR3d(ig,lg,m,iaer) = &331 k1*gIRgrid(grid_i,1,m,iaer) +&332 k2*gIRgrid(grid_i+1,1,m,iaer)709 QIRsQREF3d(ig,lg,m,iaer) = & 710 k1*qsqrefIRgrid(grid_i,1,m,iaer) + & 711 k2*qsqrefIRgrid(grid_i+1,1,m,iaer) 712 omegaIR3d(ig,lg,m,iaer) = & 713 k1*omegIRgrid(grid_i,1,m,iaer) + & 714 k2*omegIRgrid(grid_i+1,1,m,iaer) 715 gIR3d(ig,lg,m,iaer) = & 716 k1*gIRgrid(grid_i,1,m,iaer) + & 717 k2*gIRgrid(grid_i+1,1,m,iaer) 333 718 ENDDO !L_NSPECTI 334 QREFir3d(ig,lg,iaer) = & 335 k1*eprefIRgrid(grid_i,1,iaer) + & 336 k2*eprefIRgrid(grid_i+1,1,iaer) 719 QREFir3d(ig,lg,iaer) = & 720 k1*qrefIRgrid(grid_i,1,iaer) + & 721 k2*qrefIRgrid(grid_i+1,1,iaer) 722 omegaREFir3d(ig,lg,iaer) = & 723 k1*omegrefIRgrid(grid_i,1,iaer) + & 724 k2*omegrefIRgrid(grid_i+1,1,iaer) 337 725 ENDIF ! -------------------------------- 338 726 ENDDO !nlayermx 339 727 ENDDO !ngridmx 340 728 729 !================================================================== 730 731 732 341 733 ENDDO ! idomain 342 734 343 344 ENDIF ! nsize = 1 735 ENDIF ! nsize = 1 736 345 737 ENDDO ! iaer (loop on aerosol kind) 346 738 347 ! open(116,file='QIRsQREF3dO.dat')348 ! write(116,*) QIRsQREF3d349 ! close(116)350 ! open(117,file='omegaIR3dO.dat')351 ! write(117,*) omegaIR3d352 ! close(117)353 ! open(118,file='gIR3dO.dat')354 ! write(118,*) gIR3d355 ! close(118)356 ! stop357 358 739 RETURN 359 END subroutine aeroptproperties 740 END SUBROUTINE aeroptproperties 741 742 743 744 -
trunk/LMDZ.GENERIC/libf/phystd/ave_stelspec.F90
r135 r253 1 subroutine ave_stelspec( startype,STELLAR)1 subroutine ave_stelspec(STELLAR) 2 2 3 3 !================================================================== … … 28 28 29 29 #include "datafile.h" 30 #include "callkeys.h" 30 31 31 32 real*8 STELLAR(L_NSPECTV) 32 integer startype33 ! integer startype 33 34 34 35 integer Nfine … … 43 44 character(len=50) :: file_id 44 45 character(len=100) :: file_path 46 47 real lam_temp 48 double precision stel_temp 45 49 46 50 STELLAR(:)=0.0 … … 58 62 dl=lam(2)-lam(1) 59 63 60 ! load high resolution stellar data 61 if(startype.eq.1)then 62 file_id='/stellar_spectra/sol.txt' 63 tstellar=5800. 64 elseif(startype.eq.2)then 65 file_id='/stellar_spectra/gl581.txt' 66 tstellar=3200. 67 elseif(startype.eq.3)then 68 file_id='/stellar_spectra/adleo.txt' 69 tstellar=3200. 70 elseif(startype.eq.4)then 71 file_id='/stellar_spectra/gj644.txt' 72 print*,'Find out tstellar before using this star!' 73 call abort 74 elseif(startype.eq.5)then 75 file_id='/stellar_spectra/hd128167.txt' 76 tstellar=6700. ! Segura et al. (2003) 64 65 if(stelbbody)then 66 tstellar=stelTbb 67 do ifine=1,Nfine 68 lam_temp=lam(ifine) 69 call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp) 70 stel_f(ifine)=stel_temp 71 enddo 77 72 else 78 print*,'Error: unknown star type chosen' 79 call abort 73 ! load high resolution stellar data 74 if(startype.eq.1)then 75 file_id='/stellar_spectra/sol.txt' 76 tstellar=5800. 77 elseif(startype.eq.2)then 78 file_id='/stellar_spectra/gl581.txt' 79 tstellar=3200. 80 elseif(startype.eq.3)then 81 file_id='/stellar_spectra/adleo.txt' 82 tstellar=3200. 83 elseif(startype.eq.4)then 84 file_id='/stellar_spectra/gj644.txt' 85 print*,'Find out tstellar before using this star!' 86 call abort 87 elseif(startype.eq.5)then 88 file_id='/stellar_spectra/hd128167.txt' 89 tstellar=6700. ! Segura et al. (2003) 90 else 91 print*,'Error: unknown star type chosen' 92 call abort 93 endif 94 file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id)) 95 96 open(111,file=file_path,form='formatted') 97 do ifine=1,Nfine 98 read(111,*) stel_f(ifine) 99 enddo 100 close(111) 80 101 endif 81 file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))82 102 83 open(111,file=file_path,form='formatted') 84 do ifine=1,Nfine 85 read(111,*) stel_f(ifine) 86 enddo 87 close(111) 88 103 89 104 ! sum data by band 90 105 band=1 … … 104 119 STELLAR=STELLAR/sum(STELLAR) 105 120 121 106 122 end subroutine ave_stelspec -
trunk/LMDZ.GENERIC/libf/phystd/blackl.F
r135 r253 1 1 subroutine blackl(blalong,blat,blae) 2 c....................................................................... 2 3 3 implicit double precision (a-h,o-z) 4 c....................................................................... 5 c.....donnees de fond 4 5 ! physical constants 6 6 sigma=5.6693d-08 7 7 pi=datan(1.d0)*4.d0 … … 13 13 c1=h*(c**2) 14 14 c2=h*c/cbol 15 c.....fin de donnees de fond 16 c....................................................................... 15 16 17 17 blae=2.d0*pi*c1/blalong**5/(dexp(c2/blalong/blat)-1.d0) 18 c....................................................................... 18 19 19 20 return 20 21 end -
trunk/LMDZ.GENERIC/libf/phystd/calc_cpp3d.F90
r135 r253 13 13 !================================================================== 14 14 15 ! you still need to validate this equation vs something....16 17 15 implicit none 18 16 19 17 #include "comcstfi.h" 18 #include "callkeys.h" 19 #include "cpdet.h" 20 20 21 real cp0, dB2dT221 !real cp0, dB2dT2 22 22 real cppNI ! specific heat capacity at const. pressure 23 23 real rcpNI ! R / specific heat capacity 24 24 real t 25 25 real p 26 real calmol27 real bar28 26 29 calmol = 94.9784 ! 4.18/(mco2/1000) 30 bar = 100000 27 ! dummy function until I decide what to do! 31 28 32 cp0 = (7.70+5.3e-3*t - 8.3e-7*t**2)*calmol 33 dB2dT2 = 2.69e-5*t - 0.0098 34 cppNI = cp0 - t*(p/bar)*dB2dT2 35 !cppNI = 1000*(t/460)**0.35 ! Sebastian's version 36 29 cppNI = cpp 37 30 rcpNI = R/cppNI 38 31 -
trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90
r135 r253 7 7 ! Average the Rayleigh scattering in each band, weighting the 8 8 ! average by the blackbody function at temperature tstellar. 9 ! Works for an arbitrary mix of gases (currently CO2, N2 and 10 ! H2 are possible). 9 11 ! 10 12 ! Authors 11 13 ! ------- 12 14 ! Robin Wordsworth (2010) 13 ! Benjamin Charnay (2010)14 15 ! 15 16 ! Called by … … 24 25 25 26 use radinc_h, only: L_NSPECTV 26 use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, gastype,tauray, scalep27 use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep 27 28 28 29 implicit none 29 30 30 31 #include "comcstfi.h" 32 #include "callkeys.h" 33 #include "gases.h" 31 34 32 35 real*8 wl 33 integer N,Nfine,ifine 36 integer N,Nfine,ifine,igas 34 37 parameter(Nfine=500.0) 35 38 real*8 :: P0 = 9.423D+6 ! Rayleigh scattering reference pressure in pascals. 36 39 37 logical waverage38 real*8 tau const,tauvar,tausum,tauwei,bwidth,bstart40 logical typeknown 41 real*8 tauvar,tausum,tauwei,bwidth,bstart 39 42 double precision df 40 43 41 waverage=.true. ! do we perform a blackbody weighted average by band? 44 real*8 tauconsti(ngasmx) 45 real*8 tauvari(ngasmx) 42 46 43 if(waverage)then 44 if(gastype(1).eq.'CO2')then 45 print*,'Rayleigh scattering is for a CO2 atmosphere.' 46 tauconst = (8.7/g)*1.527*scalep/P0 47 elseif(gastype(1).eq.'AIR')then ! AIR = Earth air 48 print*,'Rayleigh scattering is for an Earth-like atmosphere.' 49 tauconst = (9.81/g)*8.569E-3*scalep/(P0/93.0) 47 ! tau0/p0=tau/p (Hansen 1974) 48 ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0) * 4*delta^2/(g*mugaz*lambda^4) 49 ! Then calculate tau0 = 1.16e-20 * 4*delta^2/(g*mugaz*lambda^4 [in um]) 50 ! where delta=n-1 51 52 typeknown=.false. 53 do igas=1,ngasmx 54 !if((igas.eq.vgas).or.(gfrac(igas).lt.1.e-4))then 55 if(igas.eq.vgas)then 56 print*,'Ignoring ',gnom(igas),' in Rayleigh scattering '// & 57 'as it is variable.' 58 elseif(gfrac(igas).lt.5.e-2)then 59 print*,'Ignoring ',gnom(igas),' in Rayleigh scattering '// & 60 'as its mixing ratio is less than 0.05.' 61 ! ignore variable gas in Rayleigh calculation 62 ! ignore gases of mixing ratio < 1e-4 in Rayleigh calculation 63 tauconsti(igas) = 0.0 50 64 else 51 print*,'Rayleigh scattering not defined for gas type ',gastype(n) 52 print*,'exiting.' 53 call abort 65 if(gnom(igas).eq.'CO2')then 66 tauconsti(igas) = (8.7/g)*1.527*scalep/P0 67 elseif(gnom(igas).eq.'N2_')then 68 tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0) 69 elseif(gnom(igas).eq.'H2_')then 70 tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0) 71 ! uses n=1.000132 from Optics, Fourth Edition. 72 ! around four times more opaque than CO2 73 ! around 5.5 times more opaque than air 74 elseif(gnom(igas).eq.'He_')then 75 print*,'Helium not ready yet!' 76 else 77 print*,'Error in calc_rayleigh: Gas species not recognised!' 78 call abort 79 endif 80 81 if(gfrac(igas).eq.1.0)then 82 print*,'Rayleigh scattering is for a pure ',gnom(igas),' atmosphere.' 83 typeknown=.true. 84 endif 85 54 86 endif 87 enddo 88 89 if(.not.typeknown)then 90 print*,'Rayleigh scattering is for a mixed gas atmosphere.' 91 typeknown=.true. 55 92 endif 93 56 94 57 95 do N=1,L_NSPECTV 58 96 59 if(waverage)then 60 tausum = 0.0 61 tauwei = 0.0 62 bstart = 10000.0/BWNV(N+1) 63 bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1)) 64 do ifine=1,Nfine 65 wl=bstart+dble(ifine)*bwidth/Nfine 97 tausum = 0.0 98 tauwei = 0.0 99 bstart = 10000.0/BWNV(N+1) 100 bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1)) 101 do ifine=1,Nfine 102 wl=bstart+dble(ifine)*bwidth/Nfine 66 103 67 if(gastype(1).eq.'CO2')then 68 tauvar = (1.0+0.013/wl**2)/wl**4 69 elseif(gastype(1).eq.'AIR')then 70 tauvar = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4 71 endif 104 tauvar=0.0 105 do igas=1,ngasmx 106 !if((igas.eq.vgas).or.(gfrac(igas).lt.1.e-4))then 107 if((igas.eq.vgas).or.(gfrac(igas).lt.5.e-2))then 108 ! ignore variable gas in Rayleigh calculation 109 ! ignore gases of mixing ratio < 1e-4 in Rayleigh calculation 110 tauvari(igas) = 0.0 111 else 112 if(gnom(igas).eq.'CO2')then 113 tauvari(igas) = (1.0+0.013/wl**2)/wl**4 114 elseif(gnom(igas).eq.'N2_')then 115 tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4 116 elseif(gnom(igas).eq.'H2_')then 117 tauvari(igas) = 1.0/wl**4 ! no wvl dependence of ref. index here - improve? 118 else 119 print*,'Error in calc_rayleigh: Gas species not recognised!' 120 call abort 121 endif 122 endif 72 123 73 call blackl(dble(wl*1e-6),dble(tstellar),df) 74 df=df*bwidth/Nfine 75 tauwei=tauwei+df 76 tausum=tausum+tauvar*df 124 tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas) 125 77 126 enddo 78 TAURAY(N)=tauconst*tausum/tauwei79 else80 wl = WAVEV(N)81 if(gastype(1).eq.'CO2')then82 if(N.eq.1) print*,'Rayleigh scattering is for a CO2 atmosphere.'83 TAURAY(N) = (8.7/g)*(1.527*(1.0+0.013/wl**2)/wl**4)*scalep/P084 elseif(gastype(1).eq.'AIR')then ! AIR = Earth air85 if(N.eq.1) print*,'Rayleigh scattering is for an Earth-like atmosphere.'86 TAURAY(N) = (9.81/g)*(8.569E-3*(1.0+0.0113/wl**2+0.00013/wl**4)/wl**4)*scalep/(P0/93.0)87 else88 print*,'Rayleigh scattering not defined for gas type ',gastype(n)89 print*,'exiting.'90 call abort91 endif92 endif93 127 128 call blackl(dble(wl*1e-6),dble(tstellar),df) 129 df=df*bwidth/Nfine 130 tauwei=tauwei+df 131 tausum=tausum+tauvar*df 132 133 enddo 134 TAURAY(N)=tausum/tauwei 94 135 ! we multiply by scalep here (100) because plev, which is used in optcv, 95 136 ! is in units of mBar, so we need to convert. 137 96 138 end do 97 139 140 print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-5), & 141 ' um, tau_R = ',TAURAY(L_NSPECTV-5)*1013.25 142 print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-6), & 143 ' um, tau_R = ',TAURAY(L_NSPECTV-6)*1013.25 144 98 145 end subroutine calc_rayleigh -
trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
r135 r253 3 3 ! symbols '&' in columns 73 and 6 4 4 ! 5 COMMON/callkeys/callrad,corrk,calldifv,calladj, & 6 & co2cond,callsoil, & 7 & season,diurnal,tlocked,iradia,lwrite,calllott & 8 & ,iaervar,iddist,topdustref,callstats,calleofdump & 9 & ,enertest & 10 & ,callgasvis,callnlte,callthermos,callconduct, & 11 & calleuv,callmolvis,callmoldiff,thermochem & 12 & , solarcondate, Nmix_co2, Nmix_h2o, thermoswater & 13 & , semi & 14 & , callemis & 15 & , callg2d & 16 & , linear & 17 & , ilwd & 18 & , ilwb & 19 & , ilwn & 20 & , ncouche & 21 & , alphan & 5 COMMON/callkeys/callrad,corrk,calldifv,calladj & 6 & , co2cond,callsoil & 7 & , season,diurnal,tlocked,iradia,lwrite & 8 & , iaervar,iddist,topdustref,callstats,calleofdump & 9 & , enertest & 10 & , callgasvis & 11 & , Nmix_co2, Nmix_h2o & 12 & , dusttau & 22 13 & , nonideal & 23 14 & , meanOLR & 24 15 & , specOLR & 16 & , kastprof & 17 & , noradsurf & 18 & , Tsurfk & 19 & , Tstrat & 20 & , newtonian & 21 & , tau_relax & 22 & , testradtimes & 25 23 & , rayleigh & 26 & , ozone & 24 & , stelbbody & 25 & , stelTbb & 27 26 & , tplanet & 28 27 & , startype & … … 36 35 & , rainthreshold & 37 36 & , aerofixed & 38 & , photochem,nqchem_min 37 & , szangle & 38 & , hydrology & 39 & , sourceevol & 40 & , albedosnow & 41 & , maxicethick & 42 & , Tsaldiff & 43 & , CLFfixval & 44 & , CLFvarying & 45 & , n2mixratio & 46 & , co2supsat & 47 & , pceil 39 48 40 logical callrad,corrk,calldifv,calladj,co2cond,callsoil, & 41 & season,diurnal,tlocked,lwrite,calllott & 42 & ,callstats,calleofdump & 43 & ,callgasvis,callnlte,callthermos,callconduct, & 44 & calleuv,callmolvis,callmoldiff,thermochem,thermoswater 49 logical callrad,corrk,calldifv,calladj,co2cond,callsoil & 50 & , season,diurnal,tlocked,lwrite & 51 & , callstats,calleofdump & 52 & , callgasvis 45 53 46 54 logical enertest 47 logical callemis48 logical callg2d49 logical linear50 55 logical nonideal 51 56 logical meanOLR 52 57 logical specOLR 58 logical kastprof 59 logical newtonian 60 logical testradtimes 53 61 logical rayleigh 62 logical stelbbody 54 63 logical ozone 55 64 logical nearco2cond … … 60 69 logical water,watercond,waterrain 61 70 logical aerofixed 62 logical photochem 71 logical hydrology 72 logical sourceevol 73 logical CLFvarying 74 logical noradsurf 63 75 64 76 integer iddist 65 77 integer iaervar 66 78 integer iradia 67 integer ilwd68 integer ilwb69 integer ilwn70 integer ncouche71 integer dustbin72 79 integer startype 73 integer nqchem_min74 80 75 81 real topdustref 76 real semi77 real alphan78 real solarcondate79 82 real Nmix_co2 80 83 real Nmix_h2o 84 real dusttau 81 85 real Fat1AU 86 real stelTbb 87 real Tstrat 82 88 real tplanet 83 89 real satval 84 90 real rainthreshold 91 real szangle 92 real CLFfixval 93 real n2mixratio 94 real co2supsat 95 real pceil 96 real albedosnow 97 real maxicethick 98 real Tsaldiff 99 real Tsurfk 100 real tau_relax -
trunk/LMDZ.GENERIC/libf/phystd/callsedim.F
r135 r253 1 1 SUBROUTINE callsedim(ngrid,nlay, ptimestep, 2 2 $ pplev,zlev, pt, 3 & pq, pdqfi, pdqsed,pdqs_sed,nq )3 & pq, pdqfi, pdqsed,pdqs_sed,nq,rfall) 4 4 IMPLICIT NONE 5 5 … … 29 29 30 30 #include "fisice.h" 31 c 31 32 32 c arguments: 33 33 c ---------- … … 84 84 ! (mass (kg.m-2), thickness (m), etc. 85 85 86 87 86 do l=1,nlay 88 87 do ig=1, ngrid … … 92 91 end do 93 92 94 do iq=1,nq 95 if(radius(iq).gt.1.e-9) then ! no sedimentation for gases (defined by radius=0) 93 iq=igcm_h2o_ice 96 94 97 95 ! The value of q is updated after the other parameterisations … … 101 99 ! store locally updated tracers 102 100 zqi(ig,l)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep 103 ! if (iceparty.and.(iq.eq.igcm_h2o_ice)) then104 if (iq.eq.igcm_h2o_ice) then105 c On affecte un rayon moyen aux poussieres a chaque altitude du type :106 c r(z)=r0*exp(-z/H) avec r0=0.8 micron et H=18 km.107 c ''''''''''''''''''''''''''''''''''''''''''''''''108 rfall(ig,l)=max( rice(ig,l)*1.5,rdust(ig,l) )109 c modif FranckMM pour ameliorer cycle H2O: rfall= 20 microns110 rfall(ig,l)=min(rfall(ig,l),1.e-4)111 endif112 101 enddo 113 102 enddo ! of do l=1,nlay 103 104 114 105 115 106 !======================================================================= 116 107 ! Calculate the transport due to sedimentation for each tracer 117 108 118 if (iq.eq.igcm_h2o_ice) then 119 !if (iceparty.and.(iq.eq.igcm_h2o_ice)) then 120 call newsedim(ngrid,nlay,ngrid*nlay,ptimestep, 121 & pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq) 122 else 123 call newsedim(ngrid,nlay,1,ptimestep, 124 & pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),zqi,wq) 125 endif 109 call newsedim(ngrid,nlay,ngrid*nlay,ptimestep, 110 & pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq) 126 111 127 112 !======================================================================= … … 140 125 ENDDO 141 126 142 endif ! of if(radius(iq).gt.1.e-9) 143 enddo ! of do iq=1,nq 144 145 RETURN 146 END 127 return 128 end 147 129 -
trunk/LMDZ.GENERIC/libf/phystd/comcstfi.h
r135 r253 3 3 4 4 common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg 5 common/comcstfi/avocado ,molrad,visc5 common/comcstfi/avocado!,molrad,visc 6 6 7 7 real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg 8 real avocado ,molrad,visc8 real avocado!,molrad,visc 9 9 -
trunk/LMDZ.GENERIC/libf/phystd/comdiurn.h
r135 r253 1 COMMON/comdiurn/ldiurn, 2 s sinlon(ngridmx),coslon(ngridmx), 3 s sinlat(ngridmx),coslat(ngridmx) 4 REAL sinlon,coslon,sinlat,coslat 5 LOGICAL ldiurn 1 COMMON/comdiurn/ldiurn, & 2 & sinlon(ngridmx),coslon(ngridmx), & 3 & sinlat(ngridmx),coslat(ngridmx) 4 5 real sinlon,coslon,sinlat,coslat 6 logical ldiurn -
trunk/LMDZ.GENERIC/libf/phystd/comg1d.h
r135 r253 1 c.......................................................................2 cle COMMON pour GRADS-1D3 c(Utilise pour les sorties format Grads dans la version 1D du modele)4 c 5 con peut se dire : "on ne sauvera pas plus de 1000 variables ... hein ?"6 c 1 !....................................................................... 2 ! le COMMON pour GRADS-1D 3 ! (Utilise pour les sorties format Grads dans la version 1D du modele) 4 ! 5 ! on peut se dire : "on ne sauvera pas plus de 1000 variables ... hein ?" 6 ! 7 7 INTEGER g1d_nvarmx 8 8 PARAMETER(g1d_nvarmx=1000) 9 c 10 c* g1d_nlayer ---> nombre de couches verticales11 c* g1d_nomfich ---> nom du fichier grads12 c* g1d_unitfich ---> code du fichier grads13 c* g1d_nomctl ---> nom du fichier ctl14 c* g1d_unitctl ---> code du fichier ctl15 c* g1d_premier ---> variable logique pour dire si le fichier16 cest deja ouvert17 c* g1d_irec ---> indice de derniere ecriture18 c* g1d_nvar ---> nombre de variables deja definies a la19 cderniere ecriture20 c* g1d_nomvar ---> noms des vecteurs existants21 c* g1d_dimvar ---> taille des vecteurs22 c* g1d_titrevar ---> titres des vecteurs23 c* g1d_tmp1 ---> caractere24 c* g1d_tmp2 ---> caractere25 c 9 10 ! * g1d_nlayer ---> nombre de couches verticales 11 ! * g1d_nomfich ---> nom du fichier grads 12 ! * g1d_unitfich ---> code du fichier grads 13 ! * g1d_nomctl ---> nom du fichier ctl 14 ! * g1d_unitctl ---> code du fichier ctl 15 ! * g1d_premier ---> variable logique pour dire si le fichier 16 ! est deja ouvert 17 ! * g1d_irec ---> indice de derniere ecriture 18 ! * g1d_nvar ---> nombre de variables deja definies a la 19 ! derniere ecriture 20 ! * g1d_nomvar ---> noms des vecteurs existants 21 ! * g1d_dimvar ---> taille des vecteurs 22 ! * g1d_titrevar ---> titres des vecteurs 23 ! * g1d_tmp1 ---> caractere 24 ! * g1d_tmp2 ---> caractere 25 26 26 INTEGER g1d_nlayer 27 27 CHARACTER*100 g1d_nomfich … … 43 43 integer saveG1D 44 44 45 COMMON/COMG1DI/g1d_nlayer 46 & ,g1d_unitfich 47 & ,g1d_unitctl 48 & ,g1d_irec 49 & ,g2d_irec 50 & ,g2d_appel 51 & ,g1d_nvar 45 COMMON/COMG1DI/g1d_nlayer & 46 & ,g1d_unitfich & 47 & ,g1d_unitctl & 48 & ,g1d_irec & 49 & ,g2d_irec & 50 & ,g2d_appel & 51 & ,g1d_nvar & 52 52 & ,saveG1D 53 COMMON/COMG1DC/g1d_dimvar(0:g1d_nvarmx) 54 & ,g1d_nomfich 55 & ,g1d_nomctl 56 & ,g1d_nomvar(0:g1d_nvarmx) 57 & ,g1d_titrevar(0:g1d_nvarmx) 58 & ,g1d_tmp1 59 & ,g1d_tmp2 60 COMMON/COMG1DL/g1d_premier 53 COMMON/COMG1DC/g1d_dimvar(0:g1d_nvarmx) & 54 & ,g1d_nomfich & 55 & ,g1d_nomctl & 56 & ,g1d_nomvar(0:g1d_nvarmx) & 57 & ,g1d_titrevar(0:g1d_nvarmx) & 58 & ,g1d_tmp1 & 59 & ,g1d_tmp2 60 COMMON/COMG1DL/g1d_premier & 61 61 & ,g2d_premier 62 c 63 c.......................................................................62 ! 63 !....................................................................... -
trunk/LMDZ.GENERIC/libf/phystd/comsaison.h
r135 r253 1 c-----------------------------------------------------------------------2 cINCLUDE saison.h1 !----------------------------------------------------------------------- 2 ! INCLUDE saison.h 3 3 4 COMMON/saison/callsais,isaison,dist_s ol,declin,5 $mu0(ngridmx),fract(ngridmx)4 COMMON/saison/callsais,isaison,dist_star,declin, & 5 & mu0(ngridmx),fract(ngridmx) 6 6 7 INTEGERisaison8 LOGICALcallsais9 REAL dist_sol,declin,mu0,fract7 integer isaison 8 logical callsais 9 real dist_star,declin,mu0,fract 10 10 11 c-----------------------------------------------------------------------11 !----------------------------------------------------------------------- -
trunk/LMDZ.GENERIC/libf/phystd/condens_co2cloud.F90
r135 r253 17 17 ! 18 18 ! Inputs 19 ! ------ 19 ! ------ 20 20 ! ngrid Number of vertical columns 21 21 ! nlayer Number of layers … … 56 56 #include "callkeys.h" 57 57 #include "tracer.h" 58 #include "gases.h" 58 59 59 60 !----------------------------------------------------------------------- … … 142 143 CHARACTER(LEN=20) :: tracername ! to temporarily store text 143 144 145 integer igas 146 integer,save :: igasco2=0 147 character(len=3) :: gasname 148 144 149 real reffradmin, reffradmax 150 151 real ppco2 145 152 146 153 !----------------------------------------------------------------------- … … 154 161 !reffradmin=1.e-7 155 162 !reffradmax=5.e-4 156 reffradmin=0.5e-7157 reffradmax=0.1e-3 ! FF data158 !reffradmin=0.1e-5159 !reffradmax=0.5e-3 ! JB data163 !reffradmin=0.5e-7 164 !reffradmax=0.1e-3 ! FF data 165 reffradmin=0.1e-5 166 reffradmax=0.1e-3 ! JB data 160 167 ! improve this in the future... 161 168 162 169 ! Initialisation 163 170 IF (firstcall) THEN 164 do iq=1,nqmx 165 tracername=noms(iq) 166 if (tracername.eq."co2_ice") then 167 i_co2ice=iq 168 endif 169 enddo 170 171 172 ! find CO2 ice tracer 173 do iq=1,nqmx 174 tracername=noms(iq) 175 if (tracername.eq."co2_ice") then 176 i_co2ice=iq 177 endif 178 enddo 179 180 ! find CO2 gas 181 do igas=1,ngasmx 182 gasname=gnom(igas) 183 ! gasname=noms(igas) ! was a bug 184 if (gasname.eq."CO2") then 185 igasco2=igas 186 endif 187 enddo 188 171 189 write(*,*) "condense_co2cloud: i_co2ice=",i_co2ice 172 190 … … 179 197 print*,'In condens_co2cloud: ccond=',ccond,' latcond=',latcond 180 198 181 firstcall=.false. 199 ! Prepare special treatment if gas is not pure CO2 200 !if (addn2) then 201 ! m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) 202 ! m_noco2 = 28.02E-3 ! N2 molecular mass (kg/mol) 203 ! Compute A and B coefficient use to compute 204 ! mean molecular mass Mair defined by 205 ! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2 206 ! 1/Mair = A*q(ico2) + B 207 ! A = (1/m_co2 - 1/m_noco2) 208 ! B = 1/m_noco2 209 !endif 210 211 ! Minimum CO2 mixing ratio below which mixing occurs with layer above: 212 !qco2min =0.75 213 214 firstcall=.false. 182 215 ENDIF 183 216 zcpi=1./cpp … … 217 250 ! Atmospheric condensation 218 251 252 253 ! Compute CO2 Volume mixing ratio 254 ! ------------------------------- 255 ! if (addn2) then 256 ! DO l=1,nlayer 257 ! DO ig=1,ngrid 258 ! qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep 259 !! Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2) 260 ! mmean=1/(A*qco2 +B) 261 ! vmr_co2(ig,l) = qco2*mmean/m_co2 262 ! ENDDO 263 ! ENDDO 264 ! else 265 ! DO l=1,nlayer 266 ! DO ig=1,ngrid 267 ! vmr_co2(ig,l)=0.5 268 ! ENDDO 269 ! ENDDO 270 ! end if 271 219 272 ! Forecast the atmospheric frost temperature 'ztcond' 220 273 DO l=1,nlayer 221 274 DO ig=1,ngrid 222 call get_tcond_co2(pplay(ig,l),ztcond(ig,l),1,1) 275 ppco2=gfrac(igasco2)*pplay(ig,l) 276 call get_tcond_co2(ppco2,ztcond(ig,l)) 223 277 ENDDO 224 278 ENDDO … … 234 288 print*,'Perhaps the atmosphere is collapsing on surface...?' 235 289 endif 236 stop237 238 290 END IF 239 291 ENDDO … … 275 327 ! there should be a more elegant way of doing this... 276 328 if(reff.lt.1.e-16) reff=1.e-16 277 278 329 279 330 ! update reffrad for radiative transfer … … 290 341 call stokes & 291 342 (pplev(ig,ilay),pt(ig,ilay), & 292 reff,vstokes,rho_co2) 293 343 reff,vstokes,rho_co2) 344 345 !w(ig,ilay,i_co2ice) = 0.0 294 346 w(ig,ilay,i_co2ice) = vstokes * subptimestep * & 295 347 pplev(ig,ilay)/(r*pt(ig,ilay)) … … 369 421 ! forecast of ground temperature ztsrf and frost temperature ztcondsol 370 422 DO ig=1,ngrid 371 372 call get_tcond_co2(pp lev(ig,1),ztcondsol(ig))373 423 ppco2=gfrac(igasco2)*pplay(ig,1) 424 call get_tcond_co2(ppco2,ztcondsol(ig)) 425 374 426 ztsrf(ig) = ptsrf(ig) 375 427 … … 466 518 implicit none 467 519 468 real p, tcond 469 tcond = (-3167.8)/(log(.01*p)-23.23) ! Fanale's formula 520 #include "callkeys.h" 521 522 real p, peff, tcond 523 real, parameter :: ptriple=518000.0 524 525 peff=p!/co2supsat 526 527 if(peff.lt.ptriple)then 528 tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula 529 else 530 tcond = 684.2-92.3*log(peff)+4.32*log(peff)**2 531 ! liquid-vapour transition (based on CRC handbook 2003 data) 532 endif 470 533 return 471 534 -
trunk/LMDZ.GENERIC/libf/phystd/convadj.F
r135 r253 1 SUBROUTINE convadj(ngrid,nlay,nq,ptimestep, 2 S pplay,pplev,ppopsk, 3 $ pu,pv,ph,pq, 4 $ pdufi,pdvfi,pdhfi,pdqfi, 5 $ pduadj,pdvadj,pdhadj, 6 $ pdqadj) 7 IMPLICIT NONE 8 9 c======================================================================= 10 c 11 c ajustement convectif sec 12 c on ajoute les tendances pdhfi au profil pdh avant l'ajustement 13 c SPECIAL VERSION : if one tracer is CO2, take into account the 14 c Molecular mass variation (e.g. when CO2 condense) to trigger 15 c convection F. Forget 01/2005 16 c 17 c======================================================================= 18 19 c----------------------------------------------------------------------- 20 c declarations: 21 c ------------- 1 subroutine convadj(ngrid,nlay,nq,ptimestep, 2 & pplay,pplev,ppopsk, 3 & pu,pv,ph,pq, 4 & pdufi,pdvfi,pdhfi,pdqfi, 5 & pduadj,pdvadj,pdhadj, 6 & pdqadj) 7 8 implicit none 9 10 !================================================================== 11 ! 12 ! Purpose 13 ! ------- 14 ! Calculates dry convective adjustment. If one tracer is CO2, 15 ! we take into account the molecular mass variation (e.g. when 16 ! CO2 condenses) to trigger convection (F. Forget 01/2005) 17 ! 18 ! Authors 19 ! ------- 20 ! Original author unknown. 21 ! Modif. 2005 by F. Forget. 22 ! 23 !================================================================== 24 25 ! ------------ 26 ! Declarations 27 ! ------------ 22 28 23 29 #include "dimensions.h" … … 28 34 29 35 30 c arguments: 31 c ----------36 ! Arguments 37 ! --------- 32 38 33 39 INTEGER ngrid,nlay … … 38 44 REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay) 39 45 40 c Traceurs : 46 ! Tracers 41 47 integer nq 42 48 real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) … … 44 50 45 51 46 c local: 47 c ------52 ! Local 53 ! ----- 48 54 49 55 INTEGER ig,i,l,l1,l2,jj … … 56 62 REAL zh2(ngridmx,nlayermx), zhc(ngridmx,nlayermx) 57 63 REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc 58 59 60 c Traceurs : 64 65 ! Tracers 61 66 INTEGER iq,ico2 62 67 save ico2 … … 65 70 real m_co2, m_noco2, A , B 66 71 save A, B 67 c Temporaire (for diagnostic)68 c REAL diag_alpha(ngridmx)69 72 70 73 real mtot1, mtot2 , mm1, mm2 … … 74 77 data firstcall/.true./ 75 78 79 ! for conservation test 80 real masse,cadjncons 81 76 82 EXTERNAL SCOPY 77 c 78 c-----------------------------------------------------------------------79 c initialisation: 80 c ---------------81 c 83 84 ! -------------- 85 ! Initialisation 86 ! -------------- 87 82 88 IF (firstcall) THEN 83 89 IF(ngrid.NE.ngridmx) THEN 84 90 PRINT* 85 PRINT*,'STOP dansconvadj'91 PRINT*,'STOP in convadj' 86 92 PRINT*,'ngrid =',ngrid 87 93 PRINT*,'ngridmx =',ngridmx … … 89 95 ico2=0 90 96 if (tracer) then 91 c Prepare Special treatment if one of the traceris CO2 gas97 ! Prepare Special treatment if one of the tracers is CO2 gas 92 98 do iq=1,nqmx 93 99 if (noms(iq).eq."co2") then … … 97 103 m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) 98 104 m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) 99 cCompute A and B coefficient use to compute100 cmean molecular mass Mair defined by101 c1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2102 c1/Mair = A*q(ico2) + B105 ! Compute A and B coefficient use to compute 106 ! mean molecular mass Mair defined by 107 ! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2 108 ! 1/Mair = A*q(ico2) + B 103 109 A =(1/m_co2 - 1/m_noco2) 104 110 B=1/m_noco2 … … 108 114 firstcall=.false. 109 115 ENDIF ! of IF (firstcall) 110 111 c112 c-----------------------------------------------------------------------113 c detection des profils a modifier:114 c ---------------------------------115 c si le profil est a modifier116 c (i.e. ph(niv_sup) < ph(niv_inf) )117 c alors le tableau "vtest" est mis a .TRUE. ;118 c sinon, il reste a sa valeur initiale (.FALSE.)119 c cette operation est vectorisable120 c On en profite pour copier la valeur initiale de "ph"121 c dans le champ de travail "zh"122 123 116 124 117 DO l=1,nlay … … 140 133 end if 141 134 142 143 135 CALL scopy(ngrid*nlay,zh,1,zh2,1) 144 136 CALL scopy(ngrid*nlay,zu,1,zu2,1) … … 146 138 CALL scopy(ngrid*nlay*nq,zq,1,zq2,1) 147 139 140 ! ----------------------------- 141 ! Detection of unstable columns 142 ! ----------------------------- 143 ! If ph(above) < ph(below) we set vtest=.true. 144 148 145 DO ig=1,ngrid 149 vtest(ig)=. FALSE.150 ENDDO 151 c 146 vtest(ig)=.false. 147 ENDDO 148 152 149 if (ico2.ne.0) then 153 c Special case if one of the traceris CO2 gas150 ! Special case if one of the tracers is CO2 gas 154 151 DO l=1,nlay 155 152 DO ig=1,ngrid … … 161 158 end if 162 159 163 164 160 ! Find out which grid points are convectively unstable 165 161 DO l=2,nlay 166 162 DO ig=1,ngrid 167 IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=. TRUE.163 IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true. 168 164 ENDDO 169 165 ENDDO 170 c 166 167 ! Make a list of them 171 168 jcnt=0 172 169 DO ig=1,ngrid … … 178 175 179 176 180 c-----------------------------------------------------------------------181 c Ajustement des "jcnt" profils instables indices par "jadrs": 182 c------------------------------------------------------------183 c 177 ! --------------------------------------------------------------- 178 ! Adjustment of the "jcnt" unstable profiles indicated by "jadrs" 179 ! --------------------------------------------------------------- 180 184 181 DO jj = 1, jcnt ! loop on every convective grid point 185 c 182 186 183 i = jadrs(jj) 187 184 188 c Calcul des niveaux sigma sur cette colonne 185 ! Calculate sigma in this column 189 186 DO l=1,nlay+1 190 187 sig(l)=pplev(i,l)/pplev(i,1) … … 196 193 ENDDO 197 194 l2 = 1 198 c 199 c -- boucle de sondage vers le haut 200 c 201 cins$ Loop vers le haut sur l2 195 196 ! Test loop upwards on l2 197 202 198 DO 203 c204 199 l2 = l2 + 1 205 200 IF (l2 .GT. nlay) EXIT 206 201 IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN 207 202 208 c -- l2 est le niveau le plus haut de la colonne instable 203 ! l2 is the highest level of the unstable column 209 204 210 205 l1 = l2 - 1 … … 214 209 zhm = zh2(i, l2) 215 210 if(ico2.ne.0) zqco2m = zq2(i,l2,ico2) 216 c 217 c -- boucle de sondage vers le bas218 c Loop 211 212 ! Test loop downwards 213 219 214 DO 220 c221 215 zsm = zsm + sdsig(l) 222 216 zdsm = zdsm + dsig(l) … … 230 224 end if 231 225 232 c -- doit on etendre la colonne vers le bas?233 234 down = . FALSE.235 IF (l1 . NE. 1) THEN!-- and then236 IF (zhmc . LT. zhc(i, l1-1)) THEN237 down = . TRUE.226 ! do we have to extend the column downwards? 227 228 down = .false. 229 IF (l1 .ne. 1) then !-- and then 230 IF (zhmc .lt. zhc(i, l1-1)) then 231 down = .true. 238 232 END IF 239 233 END IF 240 234 241 IF (down) THEN 235 ! this could be a problem... 236 237 if (down) then 242 238 243 239 l1 = l1 - 1 244 240 l = l1 245 241 246 ELSE247 248 c -- peut on etendre la colonne vers le haut?249 250 IF (l2 .EQ. nlay) EXIT251 252 IF (zhc(i, l2+1) .GE. zhmc) EXIT253 c 242 else 243 244 ! can we extend the column upwards? 245 246 if (l2 .eq. nlay) exit 247 248 if (zhc(i, l2+1) .ge. zhmc) exit 249 254 250 l2 = l2 + 1 255 251 l = l2 256 c 257 END IF258 c 259 cins$ End Loop 260 ENDDO 261 c 262 c -- nouveau profil : constant (valeur moyenne) 263 c 252 253 end if 254 255 enddo 256 257 ! New constant profile (average value) 258 259 264 260 zalpha=0. 265 261 zum=0. … … 276 272 endif 277 273 zh2(i, l) = zhm 278 zum=zum+dsig(l)*zu(i,l) 279 zvm=zvm+dsig(l)*zv(i,l) 274 ! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 275 zum=zum+dsig(l)*zu2(i,l) 276 zvm=zvm+dsig(l)*zv2(i,l) 277 ! zum=zum+dsig(l)*zu(i,l) 278 ! zvm=zvm+dsig(l)*zv(i,l) 280 279 do iq=1,nq 281 zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq) 280 zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq) 281 ! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq) 282 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 283 284 ! to conserve tracers/ KE, we must calculate zum, zvm and zqm using 285 ! the up-to-date column values. If we do not do this, there are cases 286 ! where convection stops at one level and starts at the next where we 287 ! can break conservation of stuff (particularly tracers) significantly. 288 282 289 end do 283 290 ENDDO … … 292 299 zalpha=1. 293 300 ELSE 294 cIF(zalpha.LT.0.) STOP301 ! IF(zalpha.LT.0.) STOP 295 302 IF(zalpha.LT.1.e-4) zalpha=1.e-4 296 303 ENDIF … … 300 307 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) 301 308 do iq=1,nq 302 czq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))309 ! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq)) 303 310 zq2(i,l,iq)=zqm(iq) 304 311 end do 305 312 ENDDO 306 c diag_alpha(i)=zalpha !temporaire307 313 if (ico2.ne.0) then 308 314 DO l=l1, l2 … … 311 317 end if 312 318 313 319 314 320 l2 = l2 + 1 315 c 316 END IF ! fin traitement instabilité entre l1 et l2. 317 c On repart pour test à partir du l2 au dessus... 318 ENDDO ! End Loop sur l2 vers le haut 319 c 320 ENDDO 321 c 321 322 END IF ! End of l1 to l2 instability treatment 323 ! We now continue to test from l2 upwards 324 325 ENDDO ! End of upwards loop on l2 326 327 328 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 329 ! check conservation 330 cadjncons=0.0 331 if(water)then 332 do l = 1, nlay 333 masse = (pplev(i,l) - pplev(i,l+1))/g 334 iq = igcm_h2o_vap 335 cadjncons = cadjncons + 336 & masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep 337 end do 338 endif 339 340 if(cadjncons.lt.-1.e-6)then 341 print*,'convadj has just crashed...' 342 print*,'i = ',i 343 print*,'l1 = ',l1 344 print*,'l2 = ',l2 345 print*,'cadjncons = ',cadjncons 346 do l = 1, nlay 347 print*,'dsig = ',dsig(l) 348 end do 349 do l = 1, nlay 350 print*,'dsig = ',dsig(l) 351 end do 352 do l = 1, nlay 353 print*,'sig = ',sig(l) 354 end do 355 do l = 1, nlay 356 print*,'pplay(ig,:) = ',pplay(i,l) 357 end do 358 do l = 1, nlay+1 359 print*,'pplev(ig,:) = ',pplev(i,l) 360 end do 361 do l = 1, nlay 362 print*,'ph(ig,:) = ',ph(i,l) 363 end do 364 do l = 1, nlay 365 print*,'ph(ig,:) = ',ph(i,l) 366 end do 367 do l = 1, nlay 368 print*,'ph(ig,:) = ',ph(i,l) 369 end do 370 do l = 1, nlay 371 print*,'zh(ig,:) = ',zh(i,l) 372 end do 373 do l = 1, nlay 374 print*,'zh2(ig,:) = ',zh2(i,l) 375 end do 376 do l = 1, nlay 377 print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap) 378 end do 379 do l = 1, nlay 380 print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap) 381 end do 382 print*,'zqm(vap) = ',zqm(igcm_h2o_vap) 383 print*,'jadrs=',jadrs 384 385 call abort 386 endif 387 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 388 389 390 ENDDO 391 322 392 DO l=1,nlay 323 393 DO ig=1,ngrid … … 338 408 end if 339 409 340 c output 410 411 ! output 341 412 ! if (ngrid.eq.1) then 342 413 ! ig=1 … … 349 420 350 421 351 RETURN352 END422 return 423 end -
trunk/LMDZ.GENERIC/libf/phystd/datafile.h
r135 r253 4 4 ! Address of the directory containing tables of data needed by the UCM 5 5 6 character (len=100) :: datafile 7 ! data datafile /'/san/home/rdword/gcm/datagcm'/ ! gnome iDataPlex 8 data datafile /'/u/rwlmd/datagcm'/ ! planetary cluster 6 character (len=300) :: datafile 7 ! data datafile /'/home/rdword/gcm/datagcm_lite'/ ! laptop 8 data datafile /'/san/home/rdword/gcm/datagcm'/ ! gnome iDataPlex 9 ! data datafile /'/u/rwlmd/datagcm'/ ! planetary cluster 9 10 !----------------------------------------------------------------------- -
trunk/LMDZ.GENERIC/libf/phystd/dimphys.h
r135 r253 8 8 integer, parameter :: nlayermx = llm 9 9 ! nsoilmx : number of subterranean layers 10 integer, parameter :: nsoilmx = 18 10 !integer, parameter :: nsoilmx = 4 ! for a test 11 integer, parameter :: nsoilmx = 18 ! for z1=0.02 cm, depth = 104.8 m 11 12 !----------------------------------------------------------------------- -
trunk/LMDZ.GENERIC/libf/phystd/dsolver.F
r135 r253 93 93 IF(XK2(N) .EQ. 0.0) GO TO 28 94 94 c IF (ABS (XK2(N)/XK(2*N-1)) .LT. 1.E-30) XK2(N)=0.0 95 95 96 IF (ABS (XK2(N)/(XK(2*N-1)+1.e-20)) .LT. 1.E-30) XK2(N)=0.0 ! For debug only (with -Ktrap=fp option) 97 96 98 97 99 28 CONTINUE -
trunk/LMDZ.GENERIC/libf/phystd/fisice.h
r135 r253 1 c-----------------------------------------------------------------------2 cINCLUDE fisice.h1 !----------------------------------------------------------------------- 2 ! INCLUDE fisice.h 3 3 4 4 COMMON/fisice/dqsf,rice,nuice,rdust,zcondicea … … 10 10 ! MODIF_JBM_08W45 11 11 REAL rdust(ngridmx,nlayermx) ! Prescribed dust radius in each layer (m) 12 cVariables used to define water ice scavenging by atmos. CO2 condensing12 ! Variables used to define water ice scavenging by atmos. CO2 condensing 13 13 REAL zcondicea(ngridmx,nlayermx) 14 c-----------------------------------------------------------------------14 !----------------------------------------------------------------------- -
trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F
r135 r253 70 70 LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI 71 71 72 NT2 = TLEV(2*L+2)*10.0D0-499 73 NT = TLEV(2*L)*10.0D0-499 72 !NT2 = TLEV(2*L+2)*10.0D0-499 73 !NT = TLEV(2*L)*10.0D0-499 74 NT = int(TLEV(2*L)*10.0D0) - NTstar+1 75 NT2 = int(TLEV(2*L+2)*10.0D0) - NTstar+1 74 76 75 77 B1(L) = (PLANCKIR(NW,NT2)-PLANCKIR(NW,NT))/DTAU(L) … … 83 85 LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI 84 86 85 NT = TLEV(2*L+1)*10.0D0-499 86 NT2 = TLEV(2*L)*10.0D0-499 87 !NT = TLEV(2*L+1)*10.0D0-499 88 !NT2 = TLEV(2*L)*10.0D0-499 89 NT = int(TLEV(2*L+1)*10.0D0) - NTstar+1 90 NT2 = int(TLEV(2*L)*10.0D0) - NTstar+1 87 91 B1(L) = (PLANCKIR(NW,NT)-PLANCKIR(NW,NT2))/DTAU(L) 88 92 B0(L) = PLANCKIR(NW,NT2) -
trunk/LMDZ.GENERIC/libf/phystd/gfluxv.F
r135 r253 66 66 67 67 68 68 69 NAYER = L_NLAYRAD 69 70 TAUMAX = L_TAUMAX !Default is 35.0 … … 84 85 W0(L) = WDEL(L)*(1.0D0-CDEL(L)**2)/FACTOR 85 86 COSBAR(L) = CDEL(L)/(1.0D0+CDEL(L)) 87 86 88 DTAU(L) = DTDEL(L)*FACTOR 87 89 TAU(L+1) = TAU(L)+DTAU(L) … … 90 92 TAUCUM(K+1) = TAUCUM(K) + (TAUCUMIN(K+1)-TAUCUMIN(K))*FACTOR 91 93 END DO 92 93 94 94 95 ! Bottom layer … … 102 103 TAUCUM(2*L+1) = TAU(L+1) 103 104 104 105 105 BSURF = RSF*UBAR0*F0PI*EXP(-MIN(TAU(L+1),TAUMAX)/UBAR0) 106 106 ! new definition of BSURF … … 109 109 ! of the radiation scattered in the forward direction 110 110 111 112 111 C WE GO WITH THE QUADRATURE APPROACH HERE. THE "SQRT(3)" factors 113 112 C ARE THE UBARV TERM. 114 113 115 114 DO L=1,L_NLAYRAD 116 117 ALPHA(L)=SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) ) 118 115 116 ALPHA(L)=SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L) ) ) 119 117 120 118 C SET OF CONSTANTS DETERMINED BY DOM 121 119 120 ! Quadrature method 122 121 G1(L) = (SQRT(3.0)*0.5)*(2.0- W0(L)*(1.0+COSBAR(L))) 123 122 G2(L) = (SQRT(3.0)*W0(L)*0.5)*(1.0-COSBAR(L)) 124 123 G3(L) = 0.5*(1.0-SQRT(3.0)*COSBAR(L)*UBAR0) 125 124 126 127 c this is Quadrature 125 ! ----- some other methods... (RDW) ------ 126 127 ! Eddington method 128 ! G1(L) = 0.25*(7.0 - W0(L)*(4.0 - 3.0*COSBAR(L))) 129 ! G2(L) = -0.25*(1.0 - W0(L)*(4.0 - 3.0*COSBAR(L))) 130 ! G3(L) = 0.25*(2.0 - 3.0*COSBAR(L)*UBAR0) 131 132 ! delta-Eddington method 133 ! G1(L) = (7.0 - 3.0*g^2 - W0(L)*(4.0 + 3.0*g) + W0(L)*g^2*(4*beta0 + 3*g)) / & 134 ! (4* (1 - g^2*() )) 0.25*(7.0 - W0(L)*(4.0 - 3.0*COSBAR(L))) 135 136 ! Hybrid modified Eddington-delta function method 137 138 ! ---------------------------------------- 139 140 c So they use Quadrature 128 141 c but the scaling is Eddington? 129 142 … … 255 268 256 269 fluxdn = fluxdn+UBAR0*F0PI*EXP(-MIN(TAUCUM(1),TAUMAX)/UBAR0) 257 270 258 271 C This is for the "special" bottom layer, where we take 259 272 C DTAU instead of DTAU/2. … … 269 282 270 283 271 272 284 C THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0 273 285 C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN … … 304 316 305 317 306 307 308 318 RETURN 309 310 311 319 END -
trunk/LMDZ.GENERIC/libf/phystd/inifis.F
r135 r253 3 3 $ plat,plon,parea, 4 4 $ prad,pg,pr,pcpp) 5 ! 5 6 use radinc_h, only : naerkind 7 8 6 9 !======================================================================= 7 10 ! … … 70 73 logical :: parameter, doubleq=.false. 71 74 75 real psurf,pN2 ! added by RW for Gliese 581d N2+CO2 76 72 77 rad=prad 73 78 daysec=pdaysec … … 78 83 rcp=r/cpp 79 84 80 81 85 avocado = 6.02214179e23 ! added by RW 82 83 86 84 87 ! -------------------------------------------------------- … … 117 120 PRINT*,'--------------------------------------------' 118 121 119 120 122 write(*,*) "Run with or without tracer transport ?" 121 123 tracer=.false. ! default value … … 168 170 write(*,*) " enertest = ",enertest 169 171 170 write(*,*) "Save EOF profiles in file 'profiles' for ",171 & "Climate Database?"172 calleofdump=.false. ! default value173 call getin("calleofdump",calleofdump)174 write(*,*) " calleofdump = ",calleofdump175 176 write(*,*) "Dust scenario:"177 iaervar=3 ! default value178 call getin("iaervar",iaervar)179 write(*,*) " iaervar = ",iaervar180 181 write(*,*) "Dust vertical distribution:"182 write(*,*) "(=1 Dust opt.deph read in startfi;",183 & " =2 Viking scenario; =3 MGS scenario,",184 & " =4 Mars Year 24 from TES assimilation)"185 iddist=3 ! default value186 call getin("iddist",iddist)187 write(*,*) " iddist = ",iddist188 189 write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"190 topdustref= 90.0 ! default value191 call getin("topdustref",topdustref)192 write(*,*) " topdustref = ",topdustref193 194 195 172 write(*,*) "call radiative transfer ?" 196 173 callrad=.true. ! default value … … 203 180 write(*,*) " corrk = ",corrk 204 181 205 ! write(*,*) "call NLTE radiative schemes ?",206 ! & "(matters only if callrad=T)"207 ! callnlte=.false. ! default value208 ! call getin("callnlte",callnlte)209 ! write(*,*) " callnlte = ",callnlte210 211 182 write(*,*) "call gaseous absorption in the visible bands?", 212 183 & "(matters only if callrad=T)" … … 242 213 write(*,*) " co2cond = ",co2cond 243 214 215 write(*,*) "CO2 supersaturation level ?" 216 co2supsat=1.0 ! default value 217 call getin("co2supsat",co2supsat) 218 write(*,*) " co2supsat = ",co2supsat 219 220 write(*,*) "Radiative timescale for Newtonian cooling ?" 221 tau_relax=30. ! default value 222 call getin("tau_relax",tau_relax) 223 write(*,*) " tau_relax = ",tau_relax 224 tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds 225 244 226 write(*,*)"call thermal conduction in the soil ?" 245 227 callsoil=.true. ! default value … … 247 229 write(*,*) " callsoil = ",callsoil 248 230 249 ! write(*,*)"call Lott's gravity wave/subgrid topography ", 250 ! & "scheme ?" 251 ! calllott=.true. ! default value 252 ! call getin("calllott",calllott) 253 ! write(*,*)" calllott = ",calllott 254 255 write(*,*)"rad.transfer is computed every iradia", 231 write(*,*)"Rad transfer is computed every iradia", 256 232 & " physical timestep" 257 233 iradia=1 ! default value … … 264 240 write(*,*)" rayleigh = ",rayleigh 265 241 266 write(*,*)"Parametrized Earth-like ozone absorption ?" 267 ozone=.false. 268 call getin("ozone",ozone) 269 write(*,*) " ozone = ",ozone 242 write(*,*) "Use blackbody for stellar spectrum ?" 243 stelbbody=.false. ! default value 244 call getin("stelbbody",stelbbody) 245 write(*,*) " stelbbody = ",stelbbody 246 247 write(*,*) "Stellar blackbody temperature ?" 248 stelTbb=5800.0 ! default value 249 call getin("stelTbb",stelTbb) 250 write(*,*) " stelTbb = ",stelTbb 270 251 271 252 write(*,*)"Output mean OLR in 1D?" … … 279 260 write(*,*)" specOLR = ",specOLR 280 261 262 write(*,*)"Operate in kastprof mode?" 263 kastprof=.false. 264 call getin("kastprof",kastprof) 265 write(*,*)" kastprof = ",kastprof 266 267 write(*,*)"Surface temperature in kastprof mode?" 268 Tsurfk=273.15 269 call getin("Tsurfk",Tsurfk) 270 write(*,*)" Tsurfk = ",Tsurfk 271 272 ! Test of incompatibility: 273 ! if kastprof used, we must be in 1D 274 if (kastprof.and.(ngridmx.gt.1)) then 275 print*,'kastprof can only be used in 1D!' 276 call abort 277 endif 278 279 write(*,*)"Stratospheric temperature for kastprof mode?" 280 Tstrat=167.0 281 call getin("Tstrat",Tstrat) 282 write(*,*)" Tstrat = ",Tstrat 283 284 write(*,*)"Remove lower boundary?" 285 noradsurf=.false. 286 call getin("noradsurf",noradsurf) 287 write(*,*)" noradsurf = ",noradsurf 288 289 ! Tests of incompatibility: 290 if (noradsurf.and.callsoil) then 291 print*,'noradsurf not compatible with soil scheme!' 292 call abort 293 endif 294 if (noradsurf.and.calldifv) then 295 print*,'noradsurf not compatible with a boundary layer!' 296 call abort 297 endif 298 299 300 write(*,*)"Use Newtonian cooling for radiative transfer?" 301 newtonian=.false. 302 call getin("newtonian",newtonian) 303 write(*,*)" newtonian = ",newtonian 304 305 ! Tests of incompatibility: 306 if (newtonian.and.corrk) then 307 print*,'newtonian not compatible with correlated-k!' 308 call abort 309 endif 310 if (newtonian.and.calladj) then 311 print*,'newtonian not compatible with adjustment!' 312 call abort 313 endif 314 if (newtonian.and.calldifv) then 315 print*,'newtonian not compatible with a boundary layer!' 316 call abort 317 endif 318 319 write(*,*)"Test physics timescale in 1D?" 320 testradtimes=.false. 321 call getin("testradtimes",testradtimes) 322 write(*,*)" testradtimes = ",testradtimes 323 324 ! Test of incompatibility: 325 ! if testradtimes used, we must be in 1D 326 if (testradtimes.and.(ngridmx.gt.1)) then 327 print*,'testradtimes can only be used in 1D!' 328 call abort 329 endif 330 281 331 write(*,*)"Default planetary temperature?" 282 332 tplanet=215.0 … … 299 349 write(*,*)" nearco2cond = ",nearco2cond 300 350 351 ! 1D solar zenith angle 352 write(*,*)"Solar zenith angle in 1D?" 353 szangle=60.0 354 call getin("szangle",szangle) 355 write(*,*)" szangle = ",szangle 356 301 357 ! TRACERS: 302 358 303 write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)" 304 dustbin=0 ! default value 305 call getin("dustbin",dustbin) 306 write(*,*)" dustbin = ",dustbin 307 308 write(*,*)"Radiatively active aerosols?" 309 aerofixed=.true. ! default value 359 write(*,*)"Fixed / inactive aerosol distribution?" 360 aerofixed=.true. ! default value 310 361 call getin("aerofixed",aerofixed) 311 362 write(*,*)" aerofixed = ",aerofixed 363 364 write(*,*)"Varying H2O cloud fraction?" 365 CLFvarying=.false. ! default value 366 call getin("CLFvarying",CLFvarying) 367 write(*,*)" CLFvarying = ",CLFvarying 368 369 write(*,*)"Value of fixed H2O cloud fraction?" 370 CLFfixval=1.0 ! default value 371 call getin("CLFfixval",CLFfixval) 372 write(*,*)" CLFfixval = ",CLFfixval 312 373 313 374 write(*,*)"Number mixing ratio of CO2 ice particles:" … … 321 382 write(*,*)" Nmix_h2o = ",Nmix_h2o 322 383 384 write(*,*)"Opacity of dust (if used):" 385 dusttau=0. ! default value 386 call getin("dusttau",dusttau) 387 write(*,*)" dusttau = ",dusttau 388 389 ! Test of incompatibility: 390 ! if less than 3 aerosols, then dusttau should = 0 391 if ((naerkind.lt.3).and.dusttau.gt.0.) then 392 print*,'Need naer>2 for radiatively active dust!' 393 stop 394 endif 395 323 396 write(*,*)"Is the variable gas species radiatively active?" 324 397 varactive=.false. … … 336 409 write(*,*)" satval = ",satval 337 410 338 ! write(*,*)"Radiatively active dust ? (matters if dustbin>0)"339 ! active=.false. ! default value340 ! call getin("active",active)341 ! write(*,*)" active = ",active342 343 411 ! Test of incompatibility: 344 412 ! if no tracers, then aerofixed should be true … … 356 424 357 425 ! Test of incompatibility: 358 ! if active is used, then dustbin should be > 0359 360 ! if (active.and.(dustbin.lt.1)) then361 ! print*,'if active is used, then dustbin should > 0'362 ! stop363 ! endif364 365 ! write(*,*)"use mass and number mixing ratios to predict",366 ! & " dust size ?"367 ! doubleq=.false. ! default value368 ! call getin("doubleq",doubleq)369 ! write(*,*)" doubleq = ",doubleq370 371 ! Test of incompatibility:372 ! if doubleq is used, then dustbin should be 1373 374 ! if (doubleq.and.(dustbin.ne.1)) then375 ! print*,'if doubleq is used, then dustbin should be 1'376 ! stop377 ! endif378 379 ! write(*,*)"dust lifted by GCM surface winds ?"380 ! lifting=.false. ! default value381 ! call getin("lifting",lifting)382 ! write(*,*)" lifting = ",lifting383 384 ! Test of incompatibility:385 ! if lifting is used, then dustbin should be > 0386 387 ! if (lifting.and.(dustbin.lt.1)) then388 ! print*,'if lifting is used, then dustbin should > 0'389 ! stop390 ! endif391 392 ! write(*,*)" dust lifted by dust devils ?"393 ! callddevil=.false. !default value394 ! call getin("callddevil",callddevil)395 ! write(*,*)" callddevil = ",callddevil396 397 398 ! Test of incompatibility:399 ! if dustdevil is used, then dustbin should be > 0400 401 ! if (callddevil.and.(dustbin.lt.1)) then402 ! print*,'if dustdevil is used, then dustbin should > 0'403 ! stop404 ! endif405 406 ! write(*,*)"Dust scavenging by CO2 snowfall ?"407 ! scavenging=.false. ! default value408 ! call getin("scavenging",scavenging)409 ! write(*,*)" scavenging = ",scavenging410 411 412 ! Test of incompatibility:413 ! if scavenging is used, then dustbin should be > 0414 415 ! if (scavenging.and.(dustbin.lt.1)) then416 ! print*,'if scavenging is used, then dustbin should > 0'417 ! stop418 ! endif419 426 420 427 write(*,*) "Gravitationnal sedimentation ?" … … 423 430 write(*,*) " sedimentation = ",sedimentation 424 431 425 ! write(*,*) "includes water ice", 426 ! & "(if true, 'water' must also be .true.)" 427 ! iceparty=.false. ! default value 428 ! call getin("iceparty",iceparty) 429 ! write(*,*) " iceparty = ",iceparty 430 431 ! write(*,*) "Radiatively active transported atmospheric ", 432 ! & "water ice ?" 433 ! activice=.false. ! default value 434 ! call getin("activice",activice) 435 ! write(*,*) " activice = ",activice 436 437 438 ! Test of incompatibility: 439 ! if activice is used, then iceparty should be used too 440 441 ! if (activice.and..not.iceparty) then 442 ! print*,'if activice is used, iceparty should be used too' 443 ! stop 444 ! endif 432 433 ! Test of incompatibility: 445 434 446 435 write(*,*) "Compute water cycle ?" … … 464 453 write(*,*) " rainthreshold = ",rainthreshold 465 454 455 write(*,*) "Include surface hydrology ?" 456 hydrology=.false. ! default value 457 call getin("hydrology",hydrology) 458 write(*,*) " hydrology = ",hydrology 459 460 write(*,*) "Evolve surface water sources ?" 461 sourceevol=.false. ! default value 462 call getin("sourceevol",sourceevol) 463 write(*,*) " sourceevol = ",sourceevol 464 465 write(*,*) "Snow albedo ?" 466 albedosnow=0.5 ! default value 467 call getin("albedosnow",albedosnow) 468 write(*,*) " albedosnow = ",albedosnow 469 470 write(*,*) "Maximum ice thickness ?" 471 maxicethick=2.0 ! default value 472 call getin("maxicethick",maxicethick) 473 write(*,*) " maxicethick = ",maxicethick 474 475 write(*,*) "Freezing point of seawater ?" 476 Tsaldiff=-1.8 ! default value 477 call getin("Tsaldiff",Tsaldiff) 478 write(*,*) " Tsaldiff = ",Tsaldiff 479 466 480 ! Test of incompatibility: 467 481 ! if watercond is used, then water should be used too … … 472 486 endif 473 487 474 write(*,*) "photochemistry: include chemical species" 475 photochem=.false. ! default value 476 call getin("photochem",photochem) 477 write(*,*) " photochem = ",photochem 478 479 480 ! THERMOSPHERE 481 482 write(*,*) "call thermosphere ?" 483 callthermos=.false. ! default value 484 call getin("callthermos",callthermos) 485 write(*,*) " callthermos = ",callthermos 486 487 write(*,*) " water included without cycle ", 488 & "(only if water=.false.)" 489 thermoswater=.false. ! default value 490 call getin("thermoswater",thermoswater) 491 write(*,*) " thermoswater = ",thermoswater 492 493 write(*,*) "call thermal conduction ?", 494 & " (only if callthermos=.true.)" 495 callconduct=.false. ! default value 496 call getin("callconduct",callconduct) 497 write(*,*) " callconduct = ",callconduct 498 499 write(*,*) "call EUV heating ?", 500 & " (only if callthermos=.true.)" 501 calleuv=.false. ! default value 502 call getin("calleuv",calleuv) 503 write(*,*) " calleuv = ",calleuv 504 505 write(*,*) "call molecular viscosity ?", 506 & " (only if callthermos=.true.)" 507 callmolvis=.false. ! default value 508 call getin("callmolvis",callmolvis) 509 write(*,*) " callmolvis = ",callmolvis 510 511 write(*,*) "call molecular diffusion ?", 512 & " (only if callthermos=.true.)" 513 callmoldiff=.false. ! default value 514 call getin("callmoldiff",callmoldiff) 515 write(*,*) " callmoldiff = ",callmoldiff 516 517 write(*,*) "call thermospheric photochemistry ?", 518 & " (only if callthermos=.true.)" 519 thermochem=.false. ! default value 520 call getin("thermochem",thermochem) 521 write(*,*) " thermochem = ",thermochem 522 523 write(*,*) "date for solar flux calculation:", 524 & " (1985 < date < 2002)" 525 write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)" 526 solarcondate=1993.4 ! default value 527 call getin("solarcondate",solarcondate) 528 write(*,*) " solarcondate = ",solarcondate 529 530 531 if (.not.callthermos) then 532 if (thermoswater) then 533 print*,'if thermoswater is set, callthermos must be true' 534 stop 535 endif 536 if (callconduct) then 537 print*,'if callconduct is set, callthermos must be true' 538 stop 539 endif 540 if (calleuv) then 541 print*,'if calleuv is set, callthermos must be true' 542 stop 543 endif 544 if (callmolvis) then 545 print*,'if callmolvis is set, callthermos must be true' 546 stop 547 endif 548 if (callmoldiff) then 549 print*,'if callmoldiff is set, callthermos must be true' 550 stop 551 endif 552 if (thermochem) then 553 print*,'if thermochem is set, callthermos must be true' 554 stop 555 endif 556 endif 557 558 ! Test of incompatibility: 559 ! if photochem is used, then water should be used too 560 561 if (photochem.and..not.water) then 562 print*,'if photochem is used, water should be used too' 563 stop 564 endif 565 566 ! if callthermos is used, then thermoswater should be used too 567 ! (if water not used already) 568 569 if (callthermos .and. .not.water) then 570 if (callthermos .and. .not.thermoswater) then 571 print*,'if callthermos is used, water or thermoswater 572 & should be used too' 573 stop 574 endif 575 endif 488 mugaz=0. 489 cpp=0. 490 call su_gases 491 call calc_cpp_mugaz 576 492 577 493 PRINT*,'--------------------------------------------' … … 594 510 PRINT*,' or each ',iradia*dtphys,' seconds' 595 511 PRINT* 596 ! -------------------------------------------------------------- 597 ! Managing the Longwave radiative transfer 598 ! -------------------------------------------------------------- 599 600 ! In most cases, the run just use the following values : 601 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 602 callemis=.true. 603 ! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01 604 ilwd=10*int(daysec/dtphys) 605 ilwn=2 606 linear=.true. 607 ncouche=3 608 alphan=0.4 609 ilwb=2 610 semi=0 611 612 c$$$! BUT people working hard on the LW may want to read them in 'radia.def' 613 c$$$! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 614 c$$$ 615 c$$$ OPEN(99,file='radia.def',status='old',form='formatted' 616 c$$$ . ,iostat=ierr) 617 c$$$ IF(ierr.EQ.0) THEN 618 c$$$ write(*,*) 'inifis: Reading radia.def !!!' 619 c$$$ READ(99,fmt='(a)') ch1 620 c$$$ READ(99,*) callemis 621 c$$$ WRITE(*,8000) ch1,callemis 622 c$$$ 623 c$$$ READ(99,fmt='(a)') ch1 624 c$$$ READ(99,*) iradia 625 c$$$ WRITE(*,8001) ch1,iradia 626 c$$$ 627 c$$$ READ(99,fmt='(a)') ch1 628 c$$$ READ(99,*) ilwd 629 c$$$ WRITE(*,8001) ch1,ilwd 630 c$$$ 631 c$$$ READ(99,fmt='(a)') ch1 632 c$$$ READ(99,*) ilwn 633 c$$$ WRITE(*,8001) ch1,ilwn 634 c$$$ 635 c$$$ READ(99,fmt='(a)') ch1 636 c$$$ READ(99,*) linear 637 c$$$ WRITE(*,8000) ch1,linear 638 c$$$ 639 c$$$ READ(99,fmt='(a)') ch1 640 c$$$ READ(99,*) ncouche 641 c$$$ WRITE(*,8001) ch1,ncouche 642 c$$$ 643 c$$$ READ(99,fmt='(a)') ch1 644 c$$$ READ(99,*) alphan 645 c$$$ WRITE(*,*) ch1,alphan 646 c$$$ 647 c$$$ READ(99,fmt='(a)') ch1 648 c$$$ READ(99,*) ilwb 649 c$$$ WRITE(*,8001) ch1,ilwb 650 c$$$ 651 c$$$ 652 c$$$ READ(99,fmt='(a)') ch1 653 c$$$ READ(99,'(l1)') callg2d 654 c$$$ WRITE(*,8000) ch1,callg2d 655 c$$$ 656 c$$$ READ(99,fmt='(a)') ch1 657 c$$$ READ(99,*) semi 658 c$$$ WRITE(*,*) ch1,semi 659 c$$$ end if 660 c$$$ CLOSE(99) 512 661 513 662 514 !----------------------------------------------------------------------- … … 678 530 pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h 679 531 680 ! managing the tracers, and tests:681 ! -------------------------------682 683 if(tracer) then684 685 ! when photochem is used, nqchem_min is the rank686 ! of the first chemical species687 688 ! Ehouarn: nqchem_min is now meaningless and no longer used689 ! nqchem_min = 1690 if (photochem .or. callthermos) then691 chem = .true.692 end if693 694 if (water .or. thermoswater)then695 h2o = .true.696 else697 h2o = .false.698 endif699 700 ! TESTS701 702 print*,'inifis: TRACERS:'703 704 c$$$ if ((doubleq).and.(h2o).and.705 c$$$ $ (chem).and.(iceparty)) then706 c$$$ print*,' 2 dust tracers (doubleq)'707 c$$$ print*,' 1 water vapour tracer'708 c$$$ print*,' 1 water ice tracer'709 c$$$ print*,nqmx-4,' chemistry tracers'710 c$$$ endif711 c$$$712 c$$$ if ((doubleq).and.(h2o).and.713 c$$$ $ (chem).and..not.(iceparty)) then714 c$$$ print*,' 2 dust tracers (doubleq)'715 c$$$ print*,' 1 water vapour tracer'716 c$$$ print*,nqmx-3,' chemistry tracers'717 c$$$ endif718 c$$$719 c$$$ if ((doubleq).and.(h2o).and.720 c$$$ $ .not.(chem).and.(iceparty)) then721 c$$$ print*,' 2 dust tracers (doubleq)'722 c$$$ print*,' 1 water vapour tracer'723 c$$$ print*,' 1 water ice tracer'724 c$$$ if (nqmx.ne.4) then725 c$$$ print*,'nqmx should be 4 with these options.'726 c$$$ print*,'(or check callphys.def)'727 c$$$ stop728 c$$$ endif729 c$$$ endif730 c$$$731 c$$$ if ((doubleq).and.(h2o).and.732 c$$$ $ .not.(chem).and..not.(iceparty)) then733 c$$$ print*,' 2 dust tracers (doubleq)'734 c$$$ print*,' 1 water vapour tracer'735 c$$$ if (nqmx.ne.3) then736 c$$$ print*,'nqmx should be 3 with these options...'737 c$$$ print*,'(or check callphys.def)'738 c$$$ stop739 c$$$ endif740 c$$$ endif741 c$$$742 c$$$ if ((doubleq).and..not.(h2o)) then743 c$$$ print*,' 2 dust tracers (doubleq)'744 c$$$ if (nqmx.ne.2) then745 c$$$ print*,'nqmx should be 2 with these options...'746 c$$$ print*,'(or check callphys.def)'747 c$$$ stop748 c$$$ endif749 c$$$ endif750 751 if (.not.(doubleq).and.(h2o).and.752 ! $ (chem).and.(iceparty)) then753 $ (chem).and.(watercond)) then754 if (dustbin.gt.0) then755 print*,dustbin,' dust bins'756 endif757 print*,nqmx-2-dustbin,' chemistry tracers'758 print*,' 1 water vapour tracer'759 print*,' 1 water ice tracer'760 endif761 762 if (.not.(doubleq).and.(h2o).and.763 ! $ (chem).and..not.(iceparty)) then764 $ (chem).and..not.(watercond)) then765 if (dustbin.gt.0) then766 print*,dustbin,' dust bins'767 endif768 print*,nqmx-1-dustbin,' chemistry tracers'769 print*,' 1 water vapour tracer'770 endif771 772 if (.not.(doubleq).and.(h2o).and.773 ! $ .not.(chem).and.(iceparty)) then774 $ .not.(chem).and.(watercond)) then775 if (dustbin.gt.0) then776 print*,dustbin,' dust bins'777 endif778 print*,' 1 water vapour tracer'779 print*,' 1 water ice tracer'780 781 if (nqmx.gt.(dustbin+2)) then782 print*,'nqmx should be ',(dustbin+2),783 $ ' with these options...'784 print*,'(or check callphys.def)'785 endif786 if (nqmx.lt.(dustbin+2)) then787 print*,dustbin,' dust bins, but this should be ok I think.'788 ! stop789 endif790 endif791 792 if (.not.(doubleq).and.(h2o).and.793 ! $ .not.(chem).and..not.(iceparty)) then794 $ .not.(chem).and..not.(watercond)) then795 if (dustbin.gt.0) then796 print*,dustbin,' dust bins'797 endif798 print*,' 1 water vapour tracer'799 if (nqmx.gt.(dustbin+1)) then800 print*,'nqmx should be ',(dustbin+1),801 $ ' with these options...'802 print*,'(or check callphys.def)'803 if (nqmx.lt.(dustbin+1)) then804 stop805 endif806 endif807 endif808 809 ! if (.not.(doubleq).and..not.(h2o)) then810 ! if (dustbin.gt.0) then811 ! print*,dustbin,' dust bins'812 ! if (nqmx.ne.dustbin) then813 ! print*,'nqmx should be ',dustbin,814 ! $ ' with these options...'815 ! print*,'(or check callphys.def)'816 ! stop817 ! endif818 ! else819 print*,'dustbin=',dustbin,820 $ ': tracer should be F with these options...'821 $ ,'UNLESS you just want to move tracers around '822 ! endif823 ! endif824 825 endif ! of if (tracer)826 827 532 RETURN 828 533 END -
trunk/LMDZ.GENERIC/libf/phystd/iniorbit.F
r135 r253 1 1 SUBROUTINE iniorbit 2 $ (pap helie,pperiheli,pyear_day,pperi_day,pobliq)2 $ (papoastr,pperiastr,pyear_day,pperi_day,pobliq) 3 3 IMPLICIT NONE 4 4 … … 13 13 c Initialisation du sous programme orbite qui calcule 14 14 c a une date donnee de l'annee de duree year_day commencant 15 c a l'equinoxe de printemps et dont le peri helie se situe15 c a l'equinoxe de printemps et dont le periastre se situe 16 16 c a la date peri_day, la distance au soleil et la declinaison. 17 17 c … … 26 26 c Input: 27 27 c ------ 28 c ap helie \ aphelie et perihelie de l'orbite29 c peri heli / en millions dekilometres.28 c apoastr \ apoastron and periastron of the orbit 29 c periastr / in millions of kilometres. 30 30 c 31 31 c======================================================================= … … 41 41 c ---------- 42 42 43 REAL pap helie,pperiheli,pyear_day,pperi_day,pobliq43 REAL papoastr,pperiastr,pyear_day,pperi_day,pobliq 44 44 45 45 c Local: … … 53 53 pi=2.*asin(1.) 54 54 55 ap helie =paphelie56 peri heli=pperiheli55 apoastr =papoastr 56 periastr=pperiastr 57 57 year_day=pyear_day 58 58 obliquit=pobliq 59 59 peri_day=pperi_day 60 60 61 PRINT*,'Peri helie en Mkm ',periheli62 PRINT*,'Ap helise en Mkm ',aphelie63 PRINT*,' obliquite en degres :',obliquit61 PRINT*,'Periastron in Mkm ',periastr 62 PRINT*,'Apoastron in Mkm ',apoastr 63 PRINT*,'Obliquity in degrees :',obliquit 64 64 unitastr=149.597927 65 e_elips=(ap helie-periheli)/(periheli+aphelie)66 p_elips=0.5*(peri heli+aphelie)*(1-e_elips*e_elips)/unitastr65 e_elips=(apoastr-periastr)/(periastr+apoastr) 66 p_elips=0.5*(periastr+apoastr)*(1-e_elips*e_elips)/unitastr 67 67 68 68 print*,'e_elips',e_elips … … 97 97 98 98 timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) 99 PRINT*,' longitude solaire du perihelietimeperi = ',timeperi99 PRINT*,'Solar longitude of periastron timeperi = ',timeperi 100 100 101 101 RETURN -
trunk/LMDZ.GENERIC/libf/phystd/initracer.F
r135 r253 62 62 63 63 64 print*,'nqmx=',nqmx65 64 ! Identify tracers by their names: (and set corresponding values of mmol) 66 65 ! 0. initialize tracer indexes to zero: … … 92 91 93 92 94 print*,'Setting dustbin = 0 in initracer.F'95 dustbin=093 !print*,'Setting dustbin = 0 in initracer.F' 94 !dustbin=0 96 95 97 96 ! 1. find dust tracers … … 222 221 c$$$ else 223 222 224 print*,dustbin 225 226 if (dustbin.gt.1) then 227 print*,'ATTENTION:', 228 $ ' properties of dust need input in initracer !!!' 229 stop 230 231 else if (dustbin.eq.1) then 232 233 c This will be used for 1 dust particle size: 234 c ------------------------------------------ 235 radius(igcm_dustbin(1))=3.e-6 236 Qext(igcm_dustbin(1))=3.04 237 alpha_lift(igcm_dustbin(1))=0.0e-6 238 alpha_devil(igcm_dustbin(1))=7.65e-9 239 qextrhor(igcm_dustbin(1))=(3./4.)*Qext(igcm_dustbin(1)) 240 & /(rho_dust*radius(igcm_dustbin(1))) 241 rho_q(igcm_dustbin(1))=rho_dust 242 243 endif 244 ! end if ! (doubleq) 223 224 c$$$ if (dustbin.gt.1) then 225 c$$$ print*,'ATTENTION:', 226 c$$$ $ ' properties of dust need input in initracer !!!' 227 c$$$ stop 228 c$$$ 229 c$$$ else if (dustbin.eq.1) then 230 c$$$ 231 c$$$c This will be used for 1 dust particle size: 232 c$$$c ------------------------------------------ 233 c$$$ radius(igcm_dustbin(1))=3.e-6 234 c$$$ Qext(igcm_dustbin(1))=3.04 235 c$$$ alpha_lift(igcm_dustbin(1))=0.0e-6 236 c$$$ alpha_devil(igcm_dustbin(1))=7.65e-9 237 c$$$ qextrhor(igcm_dustbin(1))=(3./4.)*Qext(igcm_dustbin(1)) 238 c$$$ & /(rho_dust*radius(igcm_dustbin(1))) 239 c$$$ rho_q(igcm_dustbin(1))=rho_dust 240 c$$$ 241 c$$$ endif 242 c$$$! end if ! (doubleq) 245 243 246 244 c Initialization for water vapor … … 260 258 c opacity not always realistic. 261 259 260 261 ! if(ngridmx.eq.1) 262 263 264 ! to be modified for BC+ version? 262 265 do ig=1,ngridmx 263 266 if (ngridmx.ne.1) watercaptag(ig)=.false. 264 267 dryness(ig) = 1. 265 268 enddo 269 270 271 266 272 267 273 ! IF (caps) THEN -
trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90
r135 r253 1 1 subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb) 2 2 3 3 use watercommon_h, only: To, RLVTT, RCPD … … 109 109 110 110 DO k = 1, nlayermx 111 DO i = 1, ngridmx 112 local_q(i,k) = q(i,k) 113 local_t(i,k) = t(i,k) 114 rneb(i,k) = 0.0 115 d_ql(i,k) = 0.0 116 d_t(i,k) = 0.0 117 d_q(i,k) = 0.0 118 ENDDO 111 DO i = 1, ngridmx 112 local_q(i,k) = q(i,k) 113 local_t(i,k) = t(i,k) 114 rneb(i,k) = 0.0 115 d_ql(i,k) = 0.0 116 d_t(i,k) = 0.0 117 d_q(i,k) = 0.0 118 ENDDO 119 new_qb(k)=0.0 119 120 ENDDO 120 121 … … 132 133 v_p = pplay(i,k) 133 134 134 call watersat _2(v_t,v_p,v_qs(i,k))135 call watersat(v_t,v_p,v_qs(i,k)) 135 136 call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k)) 136 137 ENDDO … … 141 142 ! DO i = 1, ngridmx 142 143 ! v_t = local_t(i,k) 143 ! IF (v_t.LT. t_coup) THEN144 ! IF (v_t.LT.To) THEN 144 145 ! print*,'RHs=',q(i,k) / v_qs(i,k) 145 146 ! ELSE … … 223 224 local_q(i,k) = new_qb(k) 224 225 local_t(i,k) = cp_new_t(k) / RCPD 225 ENDDO 226 227 ! print*,'v_qs in loop=',v_qs 228 ! print*,'v_qsd in loop=',v_qsd 229 ! print*,'new_qb in loop=',new_qb 230 ! print*,'cp_delta_t in loop=',cp_delta_t 231 ENDDO 232 226 233 227 234 !--------------------------------------------------- sounding downwards … … 242 249 ! ENDIF 243 250 244 call watersat _2(v_t,v_p,v_qs(i,k))251 call watersat(v_t,v_p,v_qs(i,k)) 245 252 call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k)) 246 253 … … 332 339 ENDDO 333 340 341 ! print*,'local_q BEFORE=',local_q 342 334 343 DO k = 1, nlayermx 335 344 DO i = 1, ngridmx … … 361 370 ! print*,'d_t=',d_t 362 371 ! print*,'d_q=',d_q 372 ! print*,'local_q=',local_q 373 ! print*,'q=',q 374 ! print*,'v_qs(i,k)=',v_qs 375 ! print*,'v_qsd(i,k)=',v_qsd 376 ! print*,'cp_delta_t(k)=',cp_delta_t 377 363 378 ! print*,'d_ql=',d_ql 364 ! print*,'i_h2o=',i_h2o 379 ! print*,'delta_q=',delta_q 380 ! print*,'zq1=',zq1 381 ! print*,'zq2=',zq2 382 !! print*,'i_h2o=',i_h2o 365 383 ! print*,'i_ice=',i_ice 366 384 ! 367 385 ! print*,'IN MANABE:' 368 386 ! print*,'d_q=',d_q … … 374 392 375 393 ! Some conservation diagnostics... 376 dEtot=0.0377 dL1tot=0.0378 dL2tot=0.0379 dqtot=0.0380 masse=0.0381 DO k = 1, nlayermx382 DO i = 1, ngridmx383 384 masse = (pplev(i,k) - pplev(i,k+1))/g385 386 dEtot = dEtot + cpp*d_t(i,k)*masse387 dL1tot = dL1tot + RLVTT*d_ql(i,k)*masse388 dL2tot = dL2tot + RLVTT*d_q(i,k)*masse ! is this line necessary?389 390 dqtot = dqtot + (d_q(i,k) + d_ql(i,k))*masse391 392 ENDDO393 ENDDO394 395 print*,'In manabe energy change=',dEtot396 print*,'In manabe condense energy change 1 =',dL1tot397 print*,'In manabe condense energy change 2 =',dL2tot398 print*,'In manabe water change=',dqtot394 ! dEtot=0.0 395 ! dL1tot=0.0 396 ! dL2tot=0.0 397 ! dqtot=0.0 398 ! masse=0.0 399 ! DO k = 1, nlayermx 400 ! DO i = 1, ngridmx 401 ! 402 ! masse = (pplev(i,k) - pplev(i,k+1))/g 403 ! 404 ! dEtot = dEtot + cpp*d_t(i,k)*masse 405 ! dL1tot = dL1tot + RLVTT*d_ql(i,k)*masse 406 ! dL2tot = dL2tot + RLVTT*d_q(i,k)*masse ! is this line necessary? 407 ! 408 ! dqtot = dqtot + (d_q(i,k) + d_ql(i,k))*masse 409 ! 410 ! ENDDO 411 ! ENDDO 412 413 ! print*,'In manabe energy change=',dEtot 414 ! print*,'In manabe condense energy change 1 =',dL1tot 415 ! print*,'In manabe condense energy change 2 =',dL2tot 416 ! print*,'In manabe water change=',dqtot 399 417 400 418 RETURN 401 419 END -
trunk/LMDZ.GENERIC/libf/phystd/newsedim.F
r135 r253 12 12 !================================================================== 13 13 14 c-----------------------------------------------------------------------15 c declarations: 16 c -------------14 !----------------------------------------------------------------------- 15 ! declarations 16 ! ------------ 17 17 18 18 #include "dimensions.h" 19 19 #include "dimphys.h" 20 20 #include "comcstfi.h" 21 c 22 c arguments: 23 c ----------21 22 ! arguments 23 ! --------- 24 24 25 25 INTEGER ngrid,nlay,naersize … … 57 57 c Physical constant 58 58 c ~~~~~~~~~~~~~~~~~ 59 !REAL visc, molrad59 REAL visc, molrad 60 60 c Gas molecular viscosity (N.s.m-2) 61 !data visc/1.e-5/ ! CO261 data visc/1.e-5/ ! CO2 62 62 c Effective gas molecular radius (m) 63 !data molrad/2.2e-10/ ! CO263 data molrad/2.2e-10/ ! CO2 64 64 65 65 c local and saved variable … … 69 69 c ** un petit test de coherence 70 70 c -------------------------- 71 72 73 !print*,'temporary fixed particle rad in newsedim!!' 71 74 72 75 IF (firstcall) THEN … … 81 84 82 85 86 83 87 !======================================================================= 84 88 ! Preliminary calculations for sedimenation velocity … … 95 99 c (correction = 0.85 for irregular particles, 0.5 for disk shaped particles) 96 100 c a = a * 0.85 101 102 97 103 ENDIF 98 99 100 101 102 104 103 105 c----------------------------------------------------------------------- … … 116 118 else 117 119 i=ngrid*(l-1)+ig 118 rfall=rd(i) 120 rfall=rd(i) ! how can this be correct!!? 119 121 endif 122 123 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 124 ! TEMPORARY MODIF BY RDW !!!! 125 !rfall=5.e-6 126 if(rfall.lt.1.e-7)then 127 rfall=1.e-7 128 endif 129 !if(rfall.gt.5.e-5)then 130 ! rfall=5.e-5 131 !endif 132 133 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 134 120 135 vstokes(ig,l) = b * rho * rfall**2 * 121 136 & (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall) -
trunk/LMDZ.GENERIC/libf/phystd/optci.F90
r135 r253 3 3 TMID,PMID,TAUGSURF,QVAR) 4 4 5 6 5 use radinc_h 7 use radcommon_h, only: gasi, tlimit, wrefVAR, Cmk,tgasref,pfgasref,wnoi 6 use radcommon_h, only: gasi, tlimit, wrefVAR, Cmk,tgasref,pfgasref,wnoi,scalep 8 7 implicit none 9 8 … … 31 30 #include "comcstfi.h" 32 31 #include "callkeys.h" 32 #include "gases.h" 33 33 34 34 … … 58 58 59 59 real*8 taugsurf(L_NSPECTI,L_NGAUSS-1) 60 real*8 dco2 60 real*8 DCONT 61 double precision wn_cont, p_cont, T_cont 61 62 62 63 ! Water mixing ratio variables … … 69 70 70 71 ! variables for k in units m^-1 71 real*8 rho, dz 72 real*8 rho, dz(L_LEVELS) 73 74 integer igas 75 76 !--- Kasting's CIA ---------------------------------------- 77 !real*8, parameter :: Ci(L_NSPECTI)=[ & 78 ! 3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7, & 79 ! 5.4E-7, 1.6E-6, 0.0, & 80 ! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 81 ! 0.0, 4.0E-7, 4.0E-6, 1.4E-5, & 82 ! 1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ] 83 !real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9, & 84 ! -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, & 85 ! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 86 ! -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ] 87 !---------------------------------------------------------- 72 88 73 89 !======================================================================= … … 75 91 ! spectral interval, NW, and each Gauss point, NG. 76 92 77 DO NG=1,L_NGAUSS-1 78 do NW=1,L_NSPECTI 79 TAUGSURF(NW,NG) = 0.0D0 80 end do 81 end do 93 taugsurf(:,:) = 0.0 94 dpr(:) = 0.0 95 lkcoef(:,:) = 0.0 82 96 83 97 do K=2,L_LEVELS 84 98 DPR(k) = PLEV(K)-PLEV(K-1) 85 86 rho = PLEV(K)/(R*TMID(K))87 dz = -DPR(k)/(g*rho)88 !print*,'rho=',rho89 !print*,'dz=',dz90 91 99 U(k) = Cmk*DPR(k) ! only Cmk line in optci.F 92 ! soon to be replaced by m^-1 !!! 100 101 !--- Kasting's CIA ---------------------------------------- 102 !dz(k)=dpr(k)*189.02*TMID(K)/(0.03720*PMID(K)) 103 ! this is CO2 path length (in cm) as written by Francois 104 ! delta_z = delta_p * R_specific * T / (g * P) 105 ! But Kasting states that W is in units of _atmosphere_ cm 106 ! So we do 107 !dz(k)=dz(k)*(PMID(K)/1013.25) 108 !dz(k)=dz(k)/100.0 ! in m for SI calc 109 !---------------------------------------------------------- 110 111 ! if we have continuum opacities, we need dz 112 dz(k)=dpr(k)*R*TMID(K)/(g*PMID(K)) 93 113 94 114 call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, & 95 115 LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K)) 96 116 97 117 do LK=1,4 98 118 LKCOEF(K,LK) = LCOEF(LK) 99 119 end do 100 120 121 101 122 DO NW=1,L_NSPECTI 102 123 do iaer=1,naerkind 103 TAEROS(K,NW,IAER) = TAUAERO(K,IAER) 124 TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER) 104 125 end do 105 126 END DO 106 127 end do ! levels 107 128 108 109 129 do K=2,L_LEVELS 110 130 do nw=1,L_NSPECTI 111 131 112 DCO2 = 0.0 ! continuum absorption (no longer used) 113 114 !if(varactive)then 115 ! call water_cont(PMID(K),TMID(K),WRATIO(K),WNOI(nw),DCO2) 116 !endif 117 118 do ng=1,L_NGAUSS-1 119 132 DCONT = 0.0 ! continuum absorption 133 134 ! include H2 continuum 135 do igas=1,ngasmx 136 if(gnom(igas).eq.'H2_')then 137 138 wn_cont = dble(wnoi(nw)) 139 T_cont = dble(TMID(k)) 140 p_cont = dble(PMID(k)*scalep*gfrac(igas)) 141 142 call interpolateH2H2(wn_cont,T_cont,p_cont,DCONT,.false.) 143 144 DCONT=DCONT*dz(k) 145 146 147 endif 148 enddo 149 150 !--- Kasting's CIA ---------------------------------------- 151 !DCO2 = dz(k)*Ci(nw)*(1.2859*PMID(k)/1000.0)*(TMID(k)/300.)**Ti(nw) 152 !DCO2 = 130*Ci(nw)*(pmid(k)/1013.25)**2*(tmid(k)/300.)**Ti(nw) * dz(k) 153 ! these two have been verified to give the same results 154 !---------------------------------------------------------- 155 156 ! Water continuum currently inactive! 157 !if(varactive)then 158 ! call water_cont(PMID(K),TMID(K),WRATIO(K),WNOI(nw),DCO2) 159 !endif 160 161 do ng=1,L_NGAUSS-1 120 162 121 163 ! Now compute TAUGAS … … 125 167 ! the water data range 126 168 127 128 if (L_REFVAR.eq.1)then ! added by RW for special no variable case 169 if(L_REFVAR.eq.1)then ! added by RW for special no variable case 129 170 KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG) 130 171 KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG) … … 133 174 else 134 175 135 KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &136 (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &137 GASI(MT(K),MP(K),NVAR(K),NW,NG))138 139 KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&140 (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &141 GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))142 143 KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&144 (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &145 GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))146 147 KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*&148 (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &149 GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))176 KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* & 177 (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - & 178 GASI(MT(K),MP(K),NVAR(K),NW,NG)) 179 180 KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* & 181 (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - & 182 GASI(MT(K),MP(K)+1,NVAR(K),NW,NG)) 183 184 KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* & 185 (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - & 186 GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG)) 187 188 KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* & 189 (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - & 190 GASI(MT(K)+1,MP(K),NVAR(K),NW,NG)) 150 191 endif 151 192 152 193 ! Interpolate the gaseous k-coefficients to the requested T,P values 153 194 154 ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) + &195 ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) + & 155 196 LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4) 156 197 157 198 TAUGAS = U(k)*ANS 158 199 159 TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS 160 161 DTAUKI(K,nw,ng) = TAUGAS 200 TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT 201 !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS 202 203 DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption 204 162 205 do iaer=1,naerkind 163 DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) & 164 + DCO2 ! For Kasting CIA 165 end do 206 DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) 207 end do ! a bug was here! 166 208 167 209 end do … … 171 213 172 214 NG = L_NGAUSS 173 DTAUKI(K,nw,ng) = 0.0 215 DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption 216 174 217 do iaer=1,naerkind 175 DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) & 176 + DCO2 ! For parameterized continuum absorption 177 end do ! a bug was found here!! 218 DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) 219 end do ! a bug was here! 178 220 179 221 end do … … 275 317 276 318 !------------------------------------------------------------------------- 277 subroutine water_cont(p,T,wratio,nu,DCO 2)319 subroutine water_cont(p,T,wratio,nu,DCONT) 278 320 ! Calculates the continuum opacity for H2O 279 ! must check somewhere that it actually _is_ H2O we are using...321 ! NOT CURRENTLY USED 280 322 281 323 implicit none 282 324 283 real p, T, wratio, nu, DCO 2325 real p, T, wratio, nu, DCONT 284 326 real h1, h2 285 327 … … 287 329 h2 = 1.25e-22 + 1.67e-19*exp(-2.62e-13*nu) 288 330 289 ! DCO2 = h1*h2*p*(p*wratio)**2/(R*T) 290 291 DCO2=0.0 ! to be implemented... 331 ! DCONT = h1*h2*p*(p*wratio)**2/(R*T) 332 ! DCONT=0.0 ! to be implemented... 292 333 293 334 return 294 335 295 336 end subroutine water_cont 337 -
trunk/LMDZ.GENERIC/libf/phystd/orbite.F
r135 r253 1 SUBROUTINE orbite(pls,pdist_sol,pdecli)2 IMPLICIT NONE1 subroutine orbite(pls,pdist_star,pdecli) 2 implicit none 3 3 4 c======================================================================= 5 c 6 c Objet: 7 c ------ 8 c 9 c Distance from sun and declimation as a function of the solar 10 c longitude Ls 11 c 12 c Interface: 13 c ---------- 14 c 15 c 16 c 17 c Arguments: 18 c ---------- 19 c 20 c Input: 21 c ------ 22 c pls Ls 23 c 24 c Output: 25 c ------- 26 c pdist_sol Distance Sun-Planet in UA 27 c pdecli declinaison ( en radians ) 28 c 29 c======================================================================= 30 c----------------------------------------------------------------------- 4 !================================================================== 5 ! 6 ! Purpose 7 ! ------- 8 ! Distance from star and declination as a function of the stellar 9 ! longitude Ls 10 ! 11 ! Inputs 12 ! ------ 13 ! pls Ls 14 ! 15 ! Outputs 16 ! ------- 17 ! pdist_star Distance Star-Planet in UA 18 ! pdecli declinaison ( in radians ) 19 ! 20 !======================================================================= 21 31 22 c Declarations: 32 23 c ------------- … … 38 29 c ---------- 39 30 40 REAL pday,pdist_s ol,pdecli,pls,i31 REAL pday,pdist_star,pdecli,pls,i 41 32 42 33 c----------------------------------------------------------------------- 43 34 44 c Distance Sun-Planet35 c Star-Planet Distance 45 36 46 pdist_s ol=p_elips/(1.+e_elips*cos(pls+timeperi))37 pdist_star = p_elips/(1.+e_elips*cos(pls+timeperi)) 47 38 48 c S olar declination39 c Stellar declination 49 40 50 41 c ********************* version before 01/01/2000 ******* 51 42 52 pdecli = asin (sin(pls)*sin(obliquit*pi/180.))43 pdecli = asin (sin(pls)*sin(obliquit*pi/180.)) 53 44 54 45 c ********************* version after 01/01/2000 ******* … … 58 49 c ****************************************************** 59 50 60 61 62 51 RETURN 63 52 END -
trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F
r135 r253 1 1 SUBROUTINE phyetat0 (fichnom,tab0,Lmodif,nsoil,nq, 2 2 . day_ini,time, 3 . tsurf,tsoil,emis,q2,qsurf )3 . tsurf,tsoil,emis,q2,qsurf,cloudfrac,totcloudfrac,hice) 4 4 implicit none 5 5 c====================================================================== … … 39 39 real qsurf(ngridmx,nq) ! tracers on surface 40 40 ! real co2ice(ngridmx) ! co2 ice cover 41 real cloudfrac(ngridmx,nlayermx) 42 real hice(ngridmx), totcloudfrac(ngridmx) 43 44 41 45 42 46 !====================================================================== … … 487 491 xmax = MAXVAL(emis) 488 492 PRINT*,'Surface emissivity <emis>:', xmin, xmax 493 494 c 495 c Cloud fraction (added by BC 2010) 496 c 497 ierr = NF_INQ_VARID (nid, "cloudfrac", nvarid) 498 IF (ierr.NE.NF_NOERR) THEN 499 PRINT*, 'phyetat0: Le champ <cloudfrac> est absent' 500 cloudfrac(:,:)=0.5 501 ! CALL abort 502 ENDIF 503 #ifdef NC_DOUBLE 504 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, cloudfrac) 505 #else 506 ierr = NF_GET_VAR_REAL(nid, nvarid, cloudfrac) 507 #endif 508 IF (ierr.NE.NF_NOERR) THEN 509 PRINT*, 'phyetat0: Lecture echouee pour <cloudfrac>' 510 CALL abort 511 ENDIF 512 xmin = 1.0E+20 513 xmax = -1.0E+20 514 xmin = MINVAL(cloudfrac) 515 xmax = MAXVAL(cloudfrac) 516 PRINT*,'Cloud fraction <cloudfrac>:', xmin, xmax 517 518 519 c 520 c Total cloud fraction (added by BC 2010) 521 c 522 ierr = NF_INQ_VARID (nid, "totcloudfrac", nvarid) 523 ! ierr = NF_INQ_VARID (nid, "totalfrac", nvarid) 524 IF (ierr.NE.NF_NOERR) THEN 525 PRINT*, 'phyetat0: Le champ <totcloudfrac> est absent' 526 totcloudfrac(:)=0.5 527 ! CALL abort 528 ENDIF 529 #ifdef NC_DOUBLE 530 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, totcloudfrac) 531 #else 532 ierr = NF_GET_VAR_REAL(nid, nvarid, totcloudfrac) 533 #endif 534 IF (ierr.NE.NF_NOERR) THEN 535 PRINT*, 'phyetat0: Lecture echouee pour <totcloudfrac>' 536 CALL abort 537 ENDIF 538 xmin = 1.0E+20 539 xmax = -1.0E+20 540 xmin = MINVAL(totcloudfrac) 541 xmax = MAXVAL(totcloudfrac) 542 PRINT*,'Cloud fraction <totcloudfrac>:', xmin, xmax 543 544 545 546 547 c 548 c Height of oceanic ice (added by BC 2010) 549 c 550 ierr = NF_INQ_VARID (nid, "hice", nvarid) 551 IF (ierr.NE.NF_NOERR) THEN 552 PRINT*, 'phyetat0: Le champ <hice> est absent' 553 CALL abort 554 ENDIF 555 #ifdef NC_DOUBLE 556 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, hice) 557 #else 558 ierr = NF_GET_VAR_REAL(nid, nvarid, hice) 559 #endif 560 IF (ierr.NE.NF_NOERR) THEN 561 PRINT*, 'phyetat0: Lecture echouee pour <hice>' 562 CALL abort 563 ENDIF 564 xmin = 1.0E+20 565 xmax = -1.0E+20 566 xmin = MINVAL(hice) 567 xmax = MAXVAL(hice) 568 PRINT*,'Height of oceanic ice <hice>:', xmin, xmax 569 570 489 571 c 490 572 c pbl wind variance … … 524 606 ELSE 525 607 txt=tnom(iq) 526 if (txt.eq."h2o_vap") then608 ! if (txt.eq."h2o_vap") then 527 609 ! There is no surface tracer for h2o_vap; 528 610 ! "h2o_ice" should be loaded instead 529 txt="h2o_ice"530 write(*,*) 'phyetat0: loading surface tracer',531 & ' h2o_ice instead of h2o_vap'532 endif611 ! txt="h2o_ice" 612 ! write(*,*) 'phyetat0: loading surface tracer', 613 ! & ' h2o_ice instead of h2o_vap' 614 ! endif 533 615 ENDIF ! of IF (oldtracernames) THEN 534 616 ierr=NF_INQ_VARID(nid,txt,nvarid) -
trunk/LMDZ.GENERIC/libf/phystd/physdem1.F
r135 r253 2 2 . phystep,day_ini, 3 3 . time,tsurf,tsoil,emis,q2,qsurf, 4 . airefi,alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe) 5 6 4 . airefi,alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe, 5 . cloudfrac,totcloudfrac,hice) 7 6 8 7 use radcommon_h, only: tauvis 9 8 10 9 implicit none 11 c------------------------------------------------------------- 12 c 13 c create physics (re-)start data file "restartfi.nc" 14 c 15 c 16 c 10 11 !================================================================== 12 ! 13 ! Purpose 14 ! ------- 15 ! Create physics (re-)start data file "restartfi.nc" 16 ! 17 ! Called by 18 ! --------- 19 ! rcm1d.F90 20 ! physiq.F90 21 ! newstart.F 22 ! 23 ! Calls 24 ! ----- 25 ! none 26 ! 27 !================================================================== 28 17 29 #include "dimensions.h" 18 30 #include "paramet.h" 19 c-----------------------------------------------------------------------20 31 #include "comvert.h" 21 32 #include "comgeom2.h" … … 39 50 integer ierr,idim1,idim2,idim3,idim4,idim5,nvarid 40 51 41 c42 52 REAL phystep,time 43 53 REAL latfi(ngridmx), lonfi(ngridmx) … … 48 58 REAL tab_cntrl(length) 49 59 50 c 60 ! added by BC 61 REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx) 62 REAL totcloudfrac(ngridmx) 63 51 64 52 65 ! EXTERNAL defrun_new,iniconst,geopot,inigeom,massdair,pression 53 66 ! EXTERNAL exner_hyb , SSUM 54 c 67 55 68 #include "serre.h" 56 69 #include "clesph0.h" … … 178 191 c Informations about Mars, only for physics 179 192 tab_cntrl(14) = year_day ! length of year (sols) ~668.6 180 tab_cntrl(15) = peri heli ! min. Sun-Marsdistance (Mkm) ~206.66181 tab_cntrl(16) = ap helie ! max. SUn-Marsdistance (Mkm) ~249.22182 tab_cntrl(17) = peri_day ! date of peri helion (sols since N. spring)193 tab_cntrl(15) = periastr ! min. star-planet distance (Mkm) ~206.66 194 tab_cntrl(16) = apoastr ! max. star-planet distance (Mkm) ~249.22 195 tab_cntrl(17) = peri_day ! date of periastron (sols since N. spring) 183 196 tab_cntrl(18) = obliquit ! Obliquity of the planet (deg) ~23.98 184 197 … … 452 465 !#endif 453 466 467 468 469 470 454 471 ! Soil Thermal inertia 455 472 … … 471 488 #endif 472 489 490 491 492 473 493 ! Surface temperature 474 494 … … 527 547 ierr = NF_REDEF (nid) 528 548 #ifdef NC_DOUBLE 529 ierr = NF_DEF_VAR (nid, 549 ierr = NF_DEF_VAR (nid,"q2", NF_DOUBLE, 2, (/idim2,idim4/),nvarid) 530 550 #else 531 551 ierr = NF_DEF_VAR (nid, "q2", NF_FLOAT, 2,(/idim2,idim4/),nvarid) … … 543 563 CALL abort 544 564 ENDIF 565 566 567 568 569 570 571 572 c cloud fraction (added by BC 2010) 573 574 ierr = NF_REDEF (nid) 575 #ifdef NC_DOUBLE 576 ierr = NF_DEF_VAR (nid, "cloudfrac", NF_DOUBLE, 1, idim2,nvarid) 577 #else 578 ierr = NF_DEF_VAR (nid, "cloudfrac", NF_FLOAT, 1, idim2,nvarid) 579 #endif 580 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14, 581 . "Cloud fraction") 582 ierr = NF_ENDDEF(nid) 583 #ifdef NC_DOUBLE 584 ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cloudfrac) 585 #else 586 ierr = NF_PUT_VAR_REAL (nid,nvarid,cloudfrac) 587 #endif 588 589 c total cloud fraction (added by BC 2010) 590 591 ierr = NF_REDEF (nid) 592 #ifdef NC_DOUBLE 593 ierr = NF_DEF_VAR (nid,"totcloudfrac", NF_DOUBLE, 1, idim2,nvarid) 594 #else 595 ierr = NF_DEF_VAR (nid, "totcloudfrac", NF_FLOAT, 1, idim2,nvarid) 596 #endif 597 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14, 598 . "Total fraction") 599 ierr = NF_ENDDEF(nid) 600 #ifdef NC_DOUBLE 601 ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,totcloudfrac) 602 #else 603 ierr = NF_PUT_VAR_REAL (nid,nvarid,totcloudfrac) 604 #endif 605 606 c oceanic ice (added by BC 2010) 607 608 ierr = NF_REDEF (nid) 609 #ifdef NC_DOUBLE 610 ierr = NF_DEF_VAR (nid, "hice", NF_DOUBLE, 1, idim2,nvarid) 611 #else 612 ierr = NF_DEF_VAR (nid, "hice", NF_FLOAT, 1, idim2,nvarid) 613 #endif 614 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21, 615 . "Height of oceanic ice") 616 ierr = NF_ENDDEF(nid) 617 #ifdef NC_DOUBLE 618 ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,hice) 619 #else 620 ierr = NF_PUT_VAR_REAL (nid,nvarid,hice) 621 #endif 622 623 624 625 626 627 628 629 630 631 545 632 546 633 c tracers … … 591 678 if ((i_h2o_vap.ne.0).and.(i_h2o_ice.eq.0)) then 592 679 ! then the index of (surface) ice is i_h2o_vap 593 i_h2o_ice=i_h2o_vap 680 print*,'water vapour but no water ice, WTF?' 681 call abort 682 i_h2o_ice=i_h2o_vap 594 683 endif 595 684 endif ! of if (.not.oldtracernames) … … 601 690 ELSE 602 691 txt=tnom(iq) 603 ! Exception: there is no water vapour surface tracer 604 if (txt.eq."h2o_vap") then 605 write(*,*)"physdem1: skipping water vapour tracer" 606 if (i_h2o_ice.eq.i_h2o_vap) then 607 ! then there is no "water ice" tracer; but still 608 ! there is some water ice on the surface 609 write(*,*)" writting water ice instead" 610 txt="h2o_ice" 611 else 612 ! there is a "water ice" tracer which has been / will be 613 ! delt with in due time 614 cycle 615 endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap) 616 endif ! of if (txt.eq."h2o_vap") 692 693 694 ! in new version, h2o_vap acts as liquid surface tracer 695 ! so the section below is now redundant 696 ! ------------------------------------------------------------------ 697 ! ! Exception: there is no water vapour surface tracer 698 ! if (txt.eq."h2o_vap") then 699 ! write(*,*)"physdem1: skipping water vapour tracer" 700 ! if (i_h2o_ice.eq.i_h2o_vap) then 701 ! ! then there is no "water ice" tracer; but still 702 ! ! there is some water ice on the surface 703 ! write(*,*)" writting water ice instead" 704 ! txt="h2o_ice" 705 ! else 706 ! ! there is a "water ice" tracer which has been / will be 707 ! ! delt with in due time 708 ! cycle 709 ! endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap) 710 ! endif ! of if (txt.eq."h2o_vap") 711 ! ------------------------------------------------------------------ 712 713 714 617 715 ENDIF ! of IF (oldtracernames) 618 716 -
trunk/LMDZ.GENERIC/libf/phystd/planete.h
r135 r253 1 c-----------------------------------------------------------------------2 cINCLUDE planet.h1 !----------------------------------------------------------------------- 2 ! INCLUDE planet.h 3 3 4 COMMON/planet/ap helie,periheli,year_day,peri_day,5 $ obliquit,nres,6 $ z0,lmixmin,emin_turb,coefvis,coefir,7 $ timeperi,e_elips,p_elips,unitastr4 COMMON/planet/apoastr,periastr,year_day,peri_day, & 5 & obliquit,nres, & 6 & z0,lmixmin,emin_turb,coefvis,coefir, & 7 & timeperi,e_elips,p_elips,unitastr 8 8 9 REAL aphelie,periheli,year_day,peri_day,10 $ obliquit,nres,11 $ z0,lmixmin,emin_turb,coefvis,coefir,12 $ timeperi,e_elips,p_elips,unitastr9 real apoastr,periastr,year_day,peri_day, & 10 & obliquit,nres, & 11 & z0,lmixmin,emin_turb,coefvis,coefir, & 12 & timeperi,e_elips,p_elips,unitastr 13 13 14 14 15 c-----------------------------------------------------------------------15 !----------------------------------------------------------------------- -
trunk/LMDZ.GENERIC/libf/phystd/profile.F
r135 r253 4 4 IMPLICIT NONE 5 5 c======================================================================= 6 c Subroutine utilisee dans le modele 1-D "testphys1d"6 c Subroutine utilisee dans "rcm1d" 7 7 c pour l'initialisation du profil atmospherique 8 8 c======================================================================= -
trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90
r135 r253 64 64 ! values manually before compiling, each time you want to change the 65 65 ! dataset. 66 66 67 67 integer, parameter :: L_NGAUSS = 17 68 68 … … 70 70 integer, parameter :: L_NTREF = 7 71 71 integer, parameter :: L_REFVAR = 7 ! earth / earlymars 72 ! CO2_H2Ovar / N2OCO2rich_H2Ovar / N2OCO2poor_H2Ovar 73 ! integer, parameter :: L_REFVAR = 4 ! N2_CH4 74 75 ! integer, parameter :: L_NPREF = 8 76 ! integer, parameter :: L_NTREF = 5 77 ! integer, parameter :: L_REFVAR = 1 ! therm_test2 78 79 ! integer, parameter :: L_NPREF = 12 80 ! integer, parameter :: L_NTREF = 20 81 ! integer, parameter :: L_REFVAR = 1 ! null (for H2 etc.) 82 83 ! integer, parameter :: L_NPREF = 10 84 ! integer, parameter :: L_NTREF = 5 85 ! integer, parameter :: L_REFVAR = 7 ! N2_CO2var 72 86 73 87 ! integer, parameter :: L_NPREF = 8 74 88 ! integer, parameter :: L_NTREF = 6 75 89 ! integer, parameter :: L_REFVAR = 1 ! pureCO2 90 91 ! integer, parameter :: L_NPREF = 7 92 ! integer, parameter :: L_NTREF = 5 93 ! integer, parameter :: L_REFVAR = 1 ! degraded_pure 76 94 77 95 ! integer, parameter :: L_NPREF = 9 … … 84 102 integer, parameter :: L_PINT = (L_NPREF-1)*5+1 85 103 integer, parameter :: NAERKIND = 2 86 integer, parameter :: L_TAUMAX = 35 104 real, parameter :: L_TAUMAX = 35 105 !integer, parameter :: L_TAUMAX = 35 106 !integer, parameter :: L_TAUMAX = 40 87 107 88 108 ! For Planck function integration: 89 109 ! equivalent temperatures are 1/10 of these values 90 110 integer, parameter :: NTstar = 500 91 integer, parameter :: NTstop = 4000 111 ! integer, parameter :: NTstop = 9000 ! for hotstuff runs 112 integer, parameter :: NTstop = 6000 ! for GJ581d / earlymars runs 113 ! integer, parameter :: NTstop = 10000 ! for H2 warming runs 92 114 93 115 ! Maximum number of grain size classes -
trunk/LMDZ.GENERIC/libf/phystd/rain.F90
r135 r253 1 2 1 subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb) 2 3 3 4 4 use watercommon_h, only: To, RLVTT, RCPD, RCPV, RV, RVTMP2 … … 39 39 REAL pplev(ngridmx,nlayermx+1) ! inter-layer pressure 40 40 REAL pplay(ngridmx,nlayermx) ! mid-layer pressure 41 REAL t(ngridmx,nlayermx) ! temperature (K) 41 REAL t(ngridmx,nlayermx) ! input temperature (K) 42 REAL zt(ngridmx,nlayermx) ! working temperature (K) 42 43 REAL ql(ngridmx,nlayermx) ! liquid water (Kg/Kg) 43 44 REAL q(ngridmx,nlayermx) ! specific humidity (Kg/Kg) … … 50 51 PARAMETER (seuil_neb=0.001) 51 52 52 REAL ct ! Inverse of cloud precipitation time53 ! REAL ct ! Inverse of cloud precipitation time 53 54 ! PARAMETER (ct=1./1800.) 54 PARAMETER (ct=1./1849.479)55 ! PARAMETER (ct=1./1849.479) 55 56 56 57 REAL cl ! Precipitation threshold … … 92 93 real tnext(ngridmx,nlayermx) 93 94 94 REAL l2c(ngridmx,nlayermx) 95 real l2c(ngridmx,nlayermx) 96 real dWtot 97 95 98 96 99 ! Indices of water vapour and water ice tracers … … 118 121 PRINT*, 'in rain.F, evap_prec=', evap_prec 119 122 120 print*,ptimestep 121 print*,1./ct 122 123 if(.not.simple)then 124 IF (ABS(ptimestep-1./ct).GT.0.001) THEN 125 PRINT*, 'Must talk to Laurent Li!!!' 126 call abort 127 ENDIF 128 endif 123 !print*,ptimestep 124 !print*,1./ct 125 !if(.not.simple)then 126 ! IF (ABS(ptimestep-1./ct).GT.0.001) THEN 127 ! PRINT*, 'Must talk to Laurent Li!!!' 128 ! call abort 129 ! ENDIF 130 !endif 129 131 130 132 firstcall = .false. … … 135 137 DO i = 1, ngridmx 136 138 137 q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap) 138 ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice) 139 zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here 140 q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep 141 ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep 142 143 !q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap) 144 !ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice) 139 145 140 146 if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water … … 167 173 DO k = 1, nlayermx 168 174 DO i = 1, ngridmx 169 ttemp = t(i,k)175 ttemp = zt(i,k) 170 176 ptemp = pplay(i,k) 171 call watersat _2(ttemp,ptemp,zqs(i,k))177 call watersat(ttemp,ptemp,zqs(i,k)) 172 178 ENDDO 173 179 ENDDO 174 180 175 ! get column / layer conversion factor (+ptimstep)181 ! get column / layer conversion factor 176 182 DO k = 1, nlayermx 177 183 DO i = 1, ngridmx 178 l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep) 184 !l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep) 185 l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g 179 186 ENDDO 180 187 ENDDO … … 189 196 DO i = 1, ngridmx 190 197 IF (zrfl(i) .GT.0.) THEN 198 199 zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! BC modif here 200 zqevt = 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & 201 *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here 202 zqevt = MAX (zqevt, 0.0) 203 zqev = MIN (zqev, zqevt) 204 zqev = MAX (zqev, 0.0) 205 zrfln(i) = zrfl(i) - zqev 206 zrfln(i) = max(zrfln(i),0.0) 207 191 208 !zqev = MAX (0.0, (zqs(i,k)-q(i,k))*eeff1 ) 192 209 !zqevt = (zrfl(i)/l2c(i,k))*eeff2 193 210 !zqev = MIN (zqev, zqevt) 194 211 !zrfln(i) = zrfl(i) - zqev*l2c(i,k) 195 196 zrfln(i) = zrfl(i) - 1.5e-5*(1.0-q(i,k)/zqs(i,k))*sqrt(zrfl(i)) 197 zrfln(i) = min(zrfln(i),zrfl(i)) 212 !zrfln(i) = zrfl(i) - 1.5e-5*(1.0-q(i,k)/zqs(i,k))*sqrt(zrfl(i)) 213 !zrfln(i) = min(zrfln(i),zrfl(i)) 198 214 ! this is what is actually written in the manual 199 215 200 d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k) 201 d_t(i,k) = d_q(i,k) * RLVTT/RCPD/(1.0+RVTMP2*q(i,k)) 202 216 d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep 217 !d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here 218 d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged! 219 203 220 zrfl(i) = zrfln(i) 204 221 ENDIF … … 214 231 215 232 DO i = 1, ngridmx 216 ttemp = t(i,k) 217 IF (ttemp .ge. To) THEN 218 lconvert=rainthreshold 219 ELSEIF (ttemp .gt. t_crit) THEN 220 lconvert=rainthreshold*(1.- t_crit/ttemp) 221 ELSE 222 lconvert=0. 223 ENDIF 224 225 IF (ql(i,k).gt.lconvert)THEN ! precipitate! 226 d_ql(i,k) = (lconvert-ql(i,k))/ptimestep 227 zrfl(i) = zrfl(i) + max(ql(i,k) - lconvert,0.0)*l2c(i,k) 228 ENDIF 233 ttemp = zt(i,k) 234 IF (ttemp .ge. To) THEN 235 lconvert=rainthreshold 236 ELSEIF (ttemp .gt. t_crit) THEN 237 lconvert=rainthreshold*(1.- t_crit/ttemp) 238 lconvert=MAX(0.0,lconvert) 239 ELSE 240 lconvert=0. 241 ENDIF 242 243 244 IF (ql(i,k).gt.1.e-9) then 245 zneb(i) = MAX(rneb(i,k), seuil_neb) 246 IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate! 247 d_ql(i,k) = -MAX((ql(i,k)-lconvert),0.0) 248 zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep 249 ENDIF 250 ENDIF 229 251 ENDDO 230 252 … … 234 256 IF (rneb(i,k).GT.0.0) THEN 235 257 zoliq(i) = ql(i,k) 236 zrho(i) = pplay(i,k) / ( t(i,k) * R )258 zrho(i) = pplay(i,k) / ( zt(i,k) * R ) 237 259 zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g) 238 zfrac(i) = ( t(i,k)-ztglace) / (To-ztglace)260 zfrac(i) = (zt(i,k)-ztglace) / (To-ztglace) 239 261 zfrac(i) = MAX(zfrac(i), 0.0) 240 262 zfrac(i) = MIN(zfrac(i), 1.0) … … 246 268 DO i = 1, ngridmx 247 269 IF (rneb(i,k).GT.0.0) THEN 248 zchau(i) = ct*ptimestep/FLOAT(ninter) * zoliq(i) &270 zchau(i) = (1./FLOAT(ninter)) * zoliq(i) & 249 271 * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) * zfrac(i) 272 ! warning: this may give dodgy results for physics calls .ne. 48 per day... 273 250 274 ! this is the ONLY place where zneb, ct and cl are used 251 275 zrhol(i) = zrho(i) * zoliq(i) / zneb(i) … … 258 282 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) 259 283 zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) 260 284 261 285 ENDIF 262 286 ENDDO … … 266 290 DO i = 1, ngridmx 267 291 IF (rneb(i,k).GT.0.0) THEN 268 d_ql(i,k) = (zoliq(i) - ql(i,k)) /ptimestep269 zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k) 292 d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep 293 zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep 270 294 ENDIF 271 295 ENDDO … … 277 301 ! Rain or snow on the ground 278 302 DO i = 1, ngridmx 303 if(zrfl(i).lt.0.0)then 304 print*,'Droplets of negative rain are falling...' 305 call abort 306 endif 279 307 IF (t(i,1) .LT. To) THEN 280 308 dqssnow(i) = zrfl(i) 309 dqsrain(i) = 0.0 281 310 ELSE 311 dqssnow(i) = 0.0 282 312 dqsrain(i) = zrfl(i) ! liquid water = ice for now 283 313 ENDIF … … 295 325 d_t(i,k) = 0.0 296 326 endif 297 dqrain(i,k,i_ice) = d_ql(i,k) 327 dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep 298 328 299 329 ENDDO 300 330 ENDDO 301 331 302 ! debugging303 ! print*,'zrfl=',zrfl304 ! print*,'dqrain=',dqrain305 ! print*,'q=',q306 ! print*,'ql=',ql307 ! print*,'dql=',d_ql308 ! DO k = 1, nlayermx309 ! DO i = 1, ngridmx310 ! if(ql(i,k).lt.0.0) then311 ! print*,'below zero!!!!'312 ! call abort313 ! endif314 ! ENDDO315 ! ENDDO316 317 318 332 RETURN 319 333 end subroutine rain -
trunk/LMDZ.GENERIC/libf/phystd/setspi.F90
r135 r253 157 157 ! force planck=sigma*eps*T^4 for each temperature in array 158 158 if(forceEC)then 159 print*,'Force F=sigma*eps*T^4 for all values of T!' 159 160 do nt=NTstar,NTstop 160 161 plancksum=0.0 -
trunk/LMDZ.GENERIC/libf/phystd/setspv.F90
r135 r253 110 110 111 111 write(*,*)'Interpolating stellar spectrum from the hires data...' 112 call ave_stelspec( startype,STELLAR)112 call ave_stelspec(STELLAR) 113 113 114 114 ! Sum the stellar flux, and write out the result. … … 120 120 write(6,'(" Stellar flux at 1 AU = ",f7.2," W m-2")') sum 121 121 print*,' ' 122 122 123 123 124 !======================================================================= -
trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F
r135 r253 32 32 real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero 33 33 34 real*8 BSURFtest ! by RW for test34 ! real*8 BSURFtest ! by RW for test 35 35 36 36 real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI) … … 45 45 46 46 NFLUXTOPI = 0.0 47 47 48 DO NW=1,L_NSPECTI 48 49 NFLUXTOPI_nu(NW) = 0.0 … … 68 69 TSURF = TLEV(L_LEVELS) 69 70 70 NTS = TSURF*10.0D0-499 71 NTT = TTOP *10.0D0-499 71 NTS = int(TSURF*10.0D0)-NTstar+1 72 NTT = int(TTOP *10.0D0)-NTstar+1 73 ! NTS = TSURF*10.0D0-499 74 ! NTT = TTOP *10.0D0-499 72 75 73 BSURFtest=0.076 ! BSURFtest=0.0 74 77 75 78 DO 501 NW=1,L_NSPECTI … … 78 81 BSURF = (1.-RSFI)*PLANCKIR(NW,NTS) ! interpolate for accuracy?? 79 82 PLTOP = PLANCKIR(NW,NTT) 80 81 !BSURFtest=BSURFtest+BSURF*DWNI(NW)82 !if(NW.eq.L_NSPECTI)then83 ! print*,'eps*sigma*T^4',5.67e-8*(1.-RSFI)*TSURF**484 ! print*,'BSURF in sfluxi=',pi*BSURFtest85 !endif86 87 83 88 84 C If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the … … 116 112 117 113 114 118 115 C NOW CALCULATE THE CUMULATIVE IR NET FLUX 119 116 … … 137 134 * (1.0-FZEROI(NW)) 138 135 139 140 136 c and same thing by spectral band... (RW) 141 137 FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + 142 138 * FMUPI(L)*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW)) 143 139 144 145 140 END DO 146 147 !fup_tmp(nw)=fup_tmp(nw)+FMUPI(L_NLEVRAD-1)*DWNI(NW)*GWEIGHT(NG)148 !fdn_tmp(nw)=fdn_tmp(nw)+FMDI(L_NLEVRAD-1)*DWNI(NW)*GWEIGHT(NG)149 !fup_tmp(nw)=fup_tmp(nw)+FMUPI(1)*DWNI(NW)*GWEIGHT(NG)150 !fdn_tmp(nw)=fdn_tmp(nw)+FMDI(1)*DWNI(NW)*GWEIGHT(NG)151 141 152 142 30 CONTINUE … … 154 144 END DO !End NGAUSS LOOP 155 145 156 ! print*,'-----------------------------------'157 !print*,'FMDI(',nw,')=',FMDI(L_NLEVRAD-1)158 !print*,'FMUPI(',nw,')=',FMUPI(L_NLEVRAD-1)159 !print*,'DWNI(',nw,')=',DWNI(nw)160 ! print*,'-----------------------------------'161 162 146 40 CONTINUE 163 147 164 148 C SPECIAL 17th Gauss point 165 166 ! print*,'fzero for ng=17',fzero167 168 149 169 150 NG = L_NGAUSS … … 203 184 END DO 204 185 205 ! print*,'-----------------------------------'206 ! print*,'nw=',nw207 ! print*,'ng=',ng208 ! print*,'FMDI=',FMDI(L_NLEVRAD-1)209 ! print*,'FMUPI=',FMUPI(L_NLEVRAD-1)210 ! print*,'-----------------------------------'211 212 186 501 CONTINUE !End Spectral Interval LOOP 213 187 214 188 C *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED**** 215 189 216 !print*,'-----------------------------------'217 !print*,'gweight=',gweight218 !print*,'FLUXDNI=',FLUXDNI(L_NLEVRAD-1)219 !print*,'FLUXUPI=',FLUXUPI(L_NLEVRAD-1)220 !print*,'-----------------------------------'221 222 !do nw=1,L_NSPECTI223 ! print*,'fup_tmp(',nw,')=',fup_tmp(nw)224 ! print*,'fdn_tmp(',nw,')=',fdn_tmp(nw)225 !enddo226 227 228 190 RETURN 229 191 END -
trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F
r135 r253 1 1 SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV, 2 * UBAR0,STEL,GWEIGHT,NFLUXTOPV,NFLUX TOPV_nu,2 * UBAR0,STEL,GWEIGHT,NFLUXTOPV,NFLUXGNDV_nu, 3 3 * FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf) 4 4 … … 19 19 real*8 NFLUXTOPV, FLUXUP, FLUXDN 20 20 real*8 NFLUXTOPV_nu(L_NSPECTV) 21 real*8 NFLUXGNDV_nu(L_NSPECTV) 21 22 real*8 GWEIGHT(L_NGAUSS) 22 23 23 24 24 integer L, NG, NW, NG1,k … … 27 27 28 28 real*8 DIFFV, DIFFVT 29 30 29 real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), fzero 31 30 … … 38 37 NFLUXTOPV = 0.0 39 38 39 DO NW=1,L_NSPECTV 40 NFLUXTOPV_nu(NW)=0.0 41 NFLUXGNDV_nu(NW)=0.0 42 END DO 40 43 41 44 DO L=1,L_NLAYRAD … … 54 57 F0PI = STEL(NW) 55 58 56 57 59 FZERO = FZEROV(NW) 58 60 IF(FZERO.ge.0.99) goto 40 … … 68 70 C SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE VISIBLE 69 71 70 BTOP = 0.0 71 72 BTOP = 0.0 73 !BSURF = 0./0. ! why was this here? 74 BSURF = 0. 72 75 C LOOP OVER THE NTERMS BEGINNING HERE 73 76 77 78 ! FACTOR = 1.0D0 - WDEL(1)*CDEL(1)**2 79 ! TAU(1) = TDEL(1)*FACTOR 80 81 74 82 ETERM = MIN(TAUV(L_NLEVRAD,NW,NG),TAUMAX) 75 83 BSURF = RSFV*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0) 76 77 84 78 85 C WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM … … 84 91 85 92 86 87 93 CALL GFLUXV(DTAUV(1,NW,NG),TAUV(1,NW,NG),TAUCUMV(1,NW,NG), 88 94 * WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV, 89 95 * BTOP,BSURF,FMUPV,FMDV,DIFFV,FLUXUP,FLUXDN) 90 91 96 92 97 C NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX … … 103 108 END DO 104 109 105 c 110 c and same thing by spectral band... (RDW) 106 111 NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW) 107 112 * +(FLUXUP-FLUXDN)*GWEIGHT(NG)* 108 113 * (1.0-FZEROV(NW)) 114 115 116 c flux at gnd (RDW) 117 NFLUXGNDV_nu(NW) = NFLUXGNDV_nu(NW) 118 * +FMDV(L_NLAYRAD)*GWEIGHT(NG)*(1.0-FZEROV(NW)) 119 109 120 110 121 C THE DIFFUSE COMPONENT OF THE DOWNWARD STELLAR FLUX … … 145 156 C NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 146 157 147 148 !print*,'fmdv',fmdv149 !print*,'fzero',fzero150 151 152 158 NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*FZERO 153 159 DO L=1,L_NLAYRAD … … 161 167 * +(FLUXUP-FLUXDN)*FZERO 162 168 169 c flux at gnd (RDW) 170 NFLUXGNDV_nu(NW) = NFLUXGNDV_nu(NW)+FMDV(L_NLAYRAD)*FZERO 171 163 172 164 173 C THE DIFFUSE COMPONENT OF THE DOWNWARD STELLAR FLUX … … 169 178 500 CONTINUE 170 179 180 171 181 C *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE VISIBLE***** 172 182 -
trunk/LMDZ.GENERIC/libf/phystd/soil.F
r135 r253 35 35 36 36 ! local saved variables: 37 ! real,save :: layer(ngridmx,nsoilmx) ! layer depth 38 real,save :: mthermdiff(ngridmx,0:nsoilmx-1) ! mid-layer thermal diffusivity 39 real,save :: thermdiff(ngridmx,nsoilmx-1) ! inter-layer thermal diffusivity 40 real,save :: coefq(0:nsoilmx-1) ! q_{k+1/2} coefficients 41 real,save :: coefd(ngridmx,nsoilmx-1) ! d_k coefficients 42 real,save :: alph(ngridmx,nsoilmx-1) ! alpha_k coefficients 43 real,save :: beta(ngridmx,nsoilmx-1) ! beta_k coefficients 44 real,save :: mu 45 37 !real,save :: mthermdiff(ngridmx,0:nsoilmx-1) ! mid-layer thermal diffusivity 38 !real,save :: thermdiff(ngridmx,nsoilmx-1) ! inter-layer thermal diffusivity 39 !real,save :: coefq(0:nsoilmx-1) ! q_{k+1/2} coefficients 40 !real,save :: coefd(ngridmx,nsoilmx-1) ! d_k coefficients 41 !real,save :: alph(ngridmx,nsoilmx-1) ! alpha_k coefficients 42 !real,save :: beta(ngridmx,nsoilmx-1) ! beta_k coefficients 43 !real,save :: mu 44 45 real mthermdiff(ngridmx,0:nsoilmx-1) ! mid-layer thermal diffusivity 46 real thermdiff(ngridmx,nsoilmx-1) ! inter-layer thermal diffusivity 47 real coefq(0:nsoilmx-1) ! q_{k+1/2} coefficients 48 real coefd(ngridmx,nsoilmx-1) ! d_k coefficients 49 real alph(ngridmx,nsoilmx-1) ! alpha_k coefficients 50 real beta(ngridmx,nsoilmx-1) ! beta_k coefficients 51 real mu 52 53 save mthermdiff,thermdiff,coefq,coefd,alph,beta,mu 54 46 55 ! local variables: 47 56 integer ig,ik 57 58 59 ! print*,'tsoil=',tsoil 60 ! print*,'tsurf=',tsurf 48 61 49 62 ! 0. Initialisations and preprocessing step … … 104 117 capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))* 105 118 & thermdiff(ig,1)/mthermdiff(ig,0)) 106 !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)119 !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig) 107 120 enddo ! of do ig=1,ngrid 108 121 109 122 else ! of if (firstcall) 123 110 124 111 125 ! 1. Compute soil temperatures … … 142 156 enddo 143 157 158 144 159 ! 3. Compute surface diffusive flux & calorific capacity 145 160 do ig=1,ngrid … … 149 164 ! & (timestep*(1.-alph(ig,1))) 150 165 ! Fstar 166 167 ! print*,'this far in soil 1' 168 ! print*,'thermdiff=',thermdiff(ig,1) 169 ! print*,'mlayer=',mlayer 170 ! print*,'beta=',beta(ig,1) 171 ! print*,'alph=',alph(ig,1) 172 ! print*,'tsoil=',tsoil(ig,1) 173 151 174 fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* 152 175 & (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1)) … … 163 186 enddo 164 187 188 165 189 end 166 190 -
trunk/LMDZ.GENERIC/libf/phystd/stokes.F90
r135 r253 30 30 31 31 ! locally saved variables 32 ! --------------------- --33 real a,b 32 ! --------------------- 33 real a,b,molrad,visc 34 34 save a,b 35 35 … … 43 43 !stop 44 44 !a = 0.707*Rgas/(4*pi*molrad**2 * avocado) 45 46 molrad=2.2e-10 ! CO2 (only used in condense_co2cloud at the moment) 47 visc=1.e-5 ! CO2 48 49 45 50 a = 0.707*R/(4*pi*molrad**2 * avocado) 46 b = (2./9.) * rho_aer * g / visc 51 b = (2./9.) * rho_aer * g / visc 52 !print*,'molrad=',molrad 53 !print*,'visc=',visc 54 !print*,'a=',a 55 !print*,'b=',b 56 !print*,'rho_aer=',rho_aer 57 !stop 58 47 59 firstcall=.false. 48 60 end if … … 52 64 ! slip-flow correction according to Rossow (Icarus 36, 1-50, 1978) 53 65 w = b * rd*rd * (1 + 1.333* (a*T/P)/rd ) 54 55 66 return 56 67 end subroutine stokes -
trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90
r135 r253 21 21 RCPD = cpp 22 22 23 RV = 1000.*R/mH2O 24 RCPV = 4. *RV 25 RVTMP2 = RCPV/RCPD-1. 23 !RV = 1000.*R/mH2O 24 RV = 1000.*8.314/mH2O ! caution! R is R/mugaz already! 26 25 26 RCPV = 4. *RV ! could be more precise... 27 !RCPV = 5/2 *8.31*1000/mH2O ! BC modif here: doesn't appear to be correct 28 29 RVTMP2 = RCPV/RCPD-1. ! not currently used... 27 30 28 31 end subroutine su_watercycle -
trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90
r135 r253 77 77 ! the single scattering parameters) 78 78 79 REAL lamref ! reference wavelengths80 REAL epref ! reference extinction ep79 REAL lamref ! reference wavelengths 80 REAL epref ! reference extinction ep 81 81 82 82 ! REAL epav(L_NSPECTI) ! Average ep (= <Qext>/Qext(lamref) if epref=1) … … 103 103 forgetit=.true. 104 104 105 ! naerkind=1 : overrides dust105 ! naerkind=1 106 106 ! visible 107 107 if(forgetit)then … … 121 121 122 122 ! NOTE: these lamref's are currently for dust, not CO2 ice. 123 ! JB thinks this shouldn't matter too much, but it needs to be fixed for the124 ! fi nal version.123 ! JB thinks this shouldn't matter too much, but it needs to be 124 ! fixed / decided for the final version. 125 125 126 126 IF (naerkind .gt. 1) THEN … … 132 132 file_id(2,2) = 'optprop_iceir_n50.dat' 133 133 lamrefir(2)=12.1E-6 ! 825cm-1 TES/MGS 134 ! ENDIF135 134 ENDIF 136 135 137 IF (naerkind .gt. 2) THEN 136 IF (naerkind .eq. 3) THEN 137 ! naerkind=3 138 ! visible 139 file_id(naerkind,1) = 'optprop_dustvis_n50.dat' 140 !lamrefvis(3)=1.5E-6 ! 1.5um OMEGA/MEx 141 lamrefvis(naerkind)=0.67e-6 142 ! infrared 143 file_id(naerkind,2) = 'optprop_dustir_n50.dat' 144 !lamrefir(3)=12.1E-6 ! 825cm-1 TES/MGS 145 lamrefir(naerkind)=9.3E-6 146 147 ENDIF 148 149 IF (naerkind .gt. 3) THEN 138 150 print*,'naerkind = ',naerkind 139 print*,'but we only have data for 2types, exiting.'151 print*,'but we only have data for 3 types, exiting.' 140 152 call abort 141 153 ENDIF … … 144 156 145 157 ! Initializations: 146 147 ! call zerophys(nsizemax*2*naerkind,radiustab)148 ! call zerophys(nsizemax*naerkind*L_NSPECTV,gvis)149 ! call zerophys(nsizemax*naerkind*L_NSPECTV,omegavis)150 ! call zerophys(nsizemax*naerkind*L_NSPECTV,QVISsQREF)151 ! call zerophys(nsizemax*naerkind*L_NSPECTI,gIR)152 ! call zerophys(nsizemax*naerkind*L_NSPECTI,omegaIR)153 ! call zerophys(nsizemax*naerkind*L_NSPECTI,QIRsQREF)154 158 155 159 radiustab(:,:,:) = 0.0 … … 308 312 endif 309 313 314 315 316 317 318 310 319 !================================================================== 311 320 ! 2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS … … 392 401 END SELECT domain 393 402 394 395 403 !======================================================================== 396 404 ! 3. Deallocate temporary variables that were read in the ASCII files … … 406 414 END DO ! Loop on idomain 407 415 !======================================================================== 416 408 417 RETURN 418 419 420 409 421 END subroutine suaer_corrk 410 422 -
trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
r135 r253 2 2 3 3 !================================================================== 4 ! 4 ! 5 5 ! Purpose 6 6 ! ------- … … 8 8 ! This subroutine is a replacement for the old 'setrad', which contained 9 9 ! both absorption and scattering data. 10 ! 10 ! 11 11 ! Authors 12 12 ! ------- 13 13 ! Adapted and generalised from the NASA Ames code by Robin Wordsworth (2009) 14 ! 14 ! 15 15 ! Summary 16 16 ! ------- 17 ! 17 ! 18 18 !================================================================== 19 19 … … 28 28 #include "datafile.h" 29 29 #include "callkeys.h" 30 #include "gases.h" 30 31 31 32 !================================================================== … … 44 45 real*8 x, xi(4), yi(4), ans, p 45 46 46 integer Nspecies 47 47 integer ngas, igas 48 49 ! temporary special for addh2 50 double precision testH2 48 51 49 52 !======================================================================= … … 51 54 file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/Q.dat' 52 55 file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id)) 56 53 57 54 58 ! check that the file exists … … 62 66 ! check that database matches varactive toggle 63 67 open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted') 64 read(111,*) Nspecies 65 66 if(Nspecies.eq.1 .and. (varactive.or.varfixed))then 68 read(111,*) ngas 69 70 if(ngas.ne.ngasmx)then 71 print*,'Number of gases in radiative transfer data (',ngas,') does not', & 72 'match that in gases.def (',ngasmx,'), exiting.' 73 call abort 74 endif 75 76 if(ngas.eq.1 .and. (varactive.or.varfixed))then 67 77 print*,'You have varactive/fixed=.true. but the database [', & 68 78 corrkdir(1:LEN_TRIM(corrkdir)), & 69 79 '] has no variable species, exiting.' 70 80 call abort 71 elseif( Nspecies.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then81 elseif(ngas.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then 72 82 print*,'You have varactive and varfixed=.false. and the database [', & 73 83 corrkdir(1:LEN_TRIM(corrkdir)), & 74 84 '] has a variable species.' 75 85 call abort 76 elseif( Nspecies.gt.3 .or. Nspecies.lt.1)then77 print*, Nspecies,' species in database [', &86 elseif(ngas.gt.4 .or. ngas.lt.1)then 87 print*,ngas,' species in database [', & 78 88 corrkdir(1:LEN_TRIM(corrkdir)), & 79 89 '], radiative code cannot handle this.' 80 90 call abort 81 endif 82 83 do n=1,Nspecies 91 endif 92 93 if(ngas.gt.3)then 94 print*,'ngas>3, EXPERIMENTAL!' 95 endif 96 97 98 do n=1,ngas 84 99 read(111,*) gastype(n) 85 100 print*,'Gas ',n,' is ',gastype(n) … … 89 104 open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted') 90 105 read(111,*) L_REFVARcheck 91 if(.not.(L_REFVARcheck.eq.L_REFVAR)) then 106 if(.not.(L_REFVARcheck.eq.L_REFVAR)) then 107 print*, L_REFVARcheck 108 print*, L_REFVAR 92 109 print*,'The size of your radiative transfer mixing ratio array does ' 93 110 print*,'not match the value given in Q.dat, exiting.' … … 96 113 read(111,*) wrefvar 97 114 close(111) 115 116 ! Check that gastype and gnom match 117 do n=1,ngas 118 print*,'Gas ',n,' is ',gnom(n) 119 if(gnom(n).ne.gastype(n))then 120 print*,'Name of a gas in radiative transfer data (',gastype(n),') does not', & 121 'match that in gases.def (',gnom(n),'), exiting.' 122 endif 123 enddo 124 print*,'Confirmed gas match in radiative transfer and gases.def!' 98 125 99 126 ! display the values 100 127 print*,'Variable gas mixing ratios:' 101 128 do n=1,L_REFVAR 102 print*,n,'.',wrefvar(n),' kg/kg' 129 !print*,n,'.',wrefvar(n),' kg/kg' ! pay attention! 130 print*,n,'.',wrefvar(n),' mol/mol' 103 131 end do 104 132 print*,'' 105 106 133 107 134 !======================================================================= … … 184 211 end do 185 212 pfgasref(L_PINT) = pgasref(L_NPREF) 186 ! Warning! this may need to be generalised if we want to use uneven grids! 187 213 print*,'Warning, pfgasref needs generalisation to uneven grids!!' 188 214 189 215 !----------------------------------------------------------------------- … … 225 251 ! Get gaseous k-coefficients and interpolate onto finer pressure grid 226 252 227 228 if (callgasvis ) then253 ! VISIBLE 254 if (callgasvis.and..not.corrkdir(1:LEN_TRIM(corrkdir)).eq.'null') then 229 255 file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat' 230 256 file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id)) … … 244 270 245 271 else 246 print*,'Visible gaseous absorption is set to zero.' 272 print*,'Visible corrk gaseous absorption is set to zero.' 273 gasv8(:,:,:,:,:)=0.0 247 274 endif 248 275 249 276 ! INFRA-RED 250 file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat' 251 file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id)) 277 if (.not.corrkdir(1:LEN_TRIM(corrkdir)).eq.'null') then 278 file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat' 279 file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id)) 252 280 253 ! check that the file exists 254 inquire(FILE=file_path,EXIST=file_ok) 255 if(.not.file_ok) then 256 write(*,*)'The file ',file_path(1:LEN_TRIM(file_path)) 257 write(*,*)'was not found by sugas_corrk.F90.' 258 write(*,*)'Are you sure you have absorption data for these bands?' 259 call abort 260 endif 281 ! check that the file exists 282 inquire(FILE=file_path,EXIST=file_ok) 283 if(.not.file_ok) then 284 write(*,*)'The file ',file_path(1:LEN_TRIM(file_path)) 285 write(*,*)'was not found by sugas_corrk.F90.' 286 write(*,*)'Are you sure you have absorption data for these bands?' 287 call abort 288 endif 289 290 open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted') 291 read(111,*) gasi8 292 close(111) 261 293 262 open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted') 263 read(111,*) gasi8 264 close(111) 265 266 do nw=1,L_NSPECTI 267 fzeroi(nw) = 0. 268 end do 269 do nw=1,L_NSPECTV 270 fzerov(nw) = 0. 271 end do 294 do nw=1,L_NSPECTI 295 fzeroi(nw) = 0. 296 end do 297 do nw=1,L_NSPECTV 298 fzerov(nw) = 0. 299 end do 300 301 else 302 print*,'Infrared corrk gaseous absorption is set to zero.' 303 gasi8(:,:,:,:,:)=0.0 304 endif 272 305 273 306 ! Take log10 of the values - this is what we will interpolate. … … 440 473 end do 441 474 475 do igas=1,ngasmx 476 if(gnom(igas).eq.'H2_')then 477 call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testH2,.true.) 478 endif 479 enddo 480 481 442 482 return 443 483 end subroutine sugas_corrk -
trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
r135 r253 102 102 c---------------------------------------------- 103 103 year_day = 669. ! length of year (sols) ~668.6 104 peri heli= 206.66 ! min. Star-Planet distance (Mkm) ~206.66105 ap helie= 249.22 ! max. Star-Planet distance (Mkm) ~249.22106 peri_day = 485. ! date of peri helion (sols since N. spring)104 periastr = 206.66 ! min. Star-Planet distance (Mkm) ~206.66 105 apoastr = 249.22 ! max. Star-Planet distance (Mkm) ~249.22 106 peri_day = 485. ! date of periastron (sols since N. spring) 107 107 obliquit = 25.19 ! Obliquity of the planet (deg) ~25.19 108 108 109 109 c additional for stokes.F added by RDW 110 110 c------------------------------------- 111 molrad=2.2e-10 ! CO2112 visc=1.e-5 ! CO2111 ! molrad=2.2e-10 ! CO2 112 ! visc=1.e-5 ! CO2 113 113 114 114 c Boundary layer and turbulence … … 180 180 c Informations about planet for the physics only 181 181 year_day = tab_cntrl(tab0+14) 182 peri heli= tab_cntrl(tab0+15)183 ap helie= tab_cntrl(tab0+16)182 periastr = tab_cntrl(tab0+15) 183 apoastr = tab_cntrl(tab0+16) 184 184 peri_day = tab_cntrl(tab0+17) 185 185 obliquit = tab_cntrl(tab0+18) … … 239 239 240 240 write(*,5) '(14) year_day',tab_cntrl(tab0+14),year_day 241 write(*,5) '(15) peri heli',tab_cntrl(tab0+15),periheli242 write(*,5) '(16) ap helie',tab_cntrl(tab0+16),aphelie241 write(*,5) '(15) periastr',tab_cntrl(tab0+15),periastr 242 write(*,5) '(16) apoastr',tab_cntrl(tab0+16),apoastr 243 243 write(*,5) '(17) peri_day',tab_cntrl(tab0+17),peri_day 244 244 write(*,5) '(18) obliquit',tab_cntrl(tab0+18),obliquit … … 288 288 write(*,*) '(27) tauvis : mean dust vis. reference ', 289 289 & 'opacity' 290 write(*,*) '(35) 291 write(*,*) '(18) 292 write(*,*) '(17) peri_day : perihelion date (solsince Ls=0)'293 write(*,*) '(15) periheli : min. sun-marsdist (Mkm)'294 write(*,*) '(16) aphelie : max. sun-mars dist (Mkm)'295 write(*,*) '(14) 290 write(*,*) '(35) volcapa : soil volumetric heat capacity' 291 write(*,*) '(18) obliquit : planet obliquity (deg)' 292 write(*,*) '(17) peri_day : periastron date (sols since Ls=0)' 293 write(*,*) '(15) periastr : min. star-planet dist (Mkm)' 294 write(*,*) '(16) apoastr : max. star-planet (Mkm)' 295 write(*,*) '(14) year_day : length of year (in sols)' 296 296 write(*,*) '(5) rad : radius of the planet (m)' 297 297 write(*,*) '(6) omeg : planet rotation rate (rad/s)' … … 445 445 write(*,*) ' peri_day (new value):',peri_day 446 446 447 else if (modif(1:len_trim(modif)) .eq. 'peri heli') then448 write(*,*) 'current value:',peri heli449 write(*,*) 'peri helion should be 206.66 on currentMars'450 write(*,*) 'enter new value:' 451 117 read(*,*,iostat=ierr) peri heli447 else if (modif(1:len_trim(modif)) .eq. 'periastr') then 448 write(*,*) 'current value:',periastr 449 write(*,*) 'periastr should be 206.66 on present-day Mars' 450 write(*,*) 'enter new value:' 451 117 read(*,*,iostat=ierr) periastr 452 452 if(ierr.ne.0) goto 117 453 453 write(*,*) 454 write(*,*) ' peri heli (new value):',periheli455 456 else if (modif(1:len_trim(modif)) .eq. 'ap helie') then457 write(*,*) 'current value:',ap helie458 write(*,*) 'ap helion should be 249.22 on currentMars'459 write(*,*) 'enter new value:' 460 118 read(*,*,iostat=ierr) ap helie454 write(*,*) ' periastr (new value):',periastr 455 456 else if (modif(1:len_trim(modif)) .eq. 'apoastr') then 457 write(*,*) 'current value:',apoastr 458 write(*,*) 'apoastr should be 249.22 on present-day Mars' 459 write(*,*) 'enter new value:' 460 118 read(*,*,iostat=ierr) apoastr 461 461 if(ierr.ne.0) goto 118 462 462 write(*,*) 463 write(*,*) ' ap helie (new value):',aphelie463 write(*,*) ' apoastr (new value):',apoastr 464 464 465 465 else if (modif(1:len_trim(modif)) .eq. 'volcapa') then … … 556 556 557 557 write(*,5) '(14) year_day',tab_cntrl(tab0+14),year_day 558 write(*,5) '(15) peri heli',tab_cntrl(tab0+15),periheli559 write(*,5) '(16) ap helie',tab_cntrl(tab0+16),aphelie558 write(*,5) '(15) periastr',tab_cntrl(tab0+15),periastr 559 write(*,5) '(16) apoastr',tab_cntrl(tab0+16),apoastr 560 560 write(*,5) '(17) peri_day',tab_cntrl(tab0+17),peri_day 561 561 write(*,5) '(18) obliquit',tab_cntrl(tab0+18),obliquit -
trunk/LMDZ.GENERIC/libf/phystd/vdifc.F
r135 r253 1 SUBROUTINE vdifc(ngrid,nlay,nq,ppopsk,2 $ ptimestep,pcapcal,lecrit,3 $pplay,pplev,pzlay,pzlev,pz0,4 $pu,pv,ph,pq,ptsrf,pemis,pqsurf,5 $pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,6 $pdudif,pdvdif,pdhdif,pdtsrf,pq2,7 $pdqdif,pdqsdif)8 9 use watercommon_h, only : RLVTT 1 subroutine vdifc(ngrid,nlay,nq,rnat,ppopsk, 2 & ptimestep,pcapcal,lecrit, 3 & pplay,pplev,pzlay,pzlev,pz0, 4 & pu,pv,ph,pq,ptsrf,pemis,pqsurf, 5 & pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf, 6 & pdudif,pdvdif,pdhdif,pdtsrf,pq2, 7 & pdqdif,pdqsdif) 8 9 use watercommon_h, only : RLVTT, To, RCPD, mx_eau_sol 10 10 11 11 implicit none 12 12 13 c======================================================================= 14 c 15 c subject: 16 c -------- 17 c Turbulent diffusion (mixing) for potential T, U, V and tracer 18 c 19 c Shema implicite 20 c On commence par rajouter au variables x la tendance physique 21 c et on resoult en fait: 22 c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1) 23 c 24 c author: 25 c ------ 26 c Hourdin/Forget/Fournier 27 c======================================================================= 28 29 c----------------------------------------------------------------------- 30 c declarations: 31 c ------------- 13 !================================================================== 14 ! 15 ! Purpose 16 ! ------- 17 ! Turbulent diffusion (mixing) for pot. T, U, V and tracers 18 ! 19 ! Implicit scheme 20 ! We start by adding to variables x the physical tendencies 21 ! already computed. We resolve the equation: 22 ! 23 ! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1) 24 ! 25 ! Authors 26 ! ------- 27 ! F. Hourdin, F. Forget, R. Fournier (199X) 28 ! R. Wordsworth, B. Charnay (2010) 29 ! 30 !================================================================== 31 32 !----------------------------------------------------------------------- 33 ! declarations 34 ! ------------ 32 35 33 36 #include "dimensions.h" … … 40 43 41 44 #include "watercap.h" 42 c 43 c arguments: 44 c ---------- 45 45 46 ! arguments 47 ! --------- 46 48 INTEGER ngrid,nlay 47 49 REAL ptimestep … … 55 57 REAL pdtsrf(ngrid),pcapcal(ngrid) 56 58 REAL pq2(ngrid,nlay+1) 57 58 c Argument added for condensation: 59 ! REAL co2ice (ngrid) 59 60 integer rnat(ngrid) 61 62 ! Arguments added for condensation 60 63 REAL ppopsk(ngrid,nlay) 61 64 logical lecrit 62 65 REAL pz0 63 66 64 c Traceurs : 67 ! Tracers 68 ! -------- 65 69 integer nq 66 REALpqsurf(ngrid,nq)70 real pqsurf(ngrid,nq) 67 71 real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 68 72 real pdqdif(ngrid,nlay,nq) 69 73 real pdqsdif(ngrid,nq) 70 74 71 c local: 72 c ------ 73 74 INTEGER ilev,ig,ilay,nlev 75 ! local 76 ! ----- 77 integer ilev,ig,ilay,nlev 75 78 76 79 REAL z4st,zdplanck(ngridmx) … … 86 89 REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx) 87 90 REAL zcst1 88 REAL zu2 89 90 EXTERNAL SSUM,SCOPY 91 REAL SSUM 91 REAL zu2!, a 92 REAL zcq(ngridmx,nlayermx),zdq(ngridmx,nlayermx) 93 REAL evap(ngridmx) 94 REAL zcq0(ngridmx),zdq0(ngridmx) 95 REAL zx_alf1(ngridmx),zx_alf2(ngridmx) 96 92 97 LOGICAL firstcall 93 98 SAVE firstcall 94 99 95 c variable added for CO2 condensation: 96 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 97 REAL hh , zhcond(ngridmx,nlayermx)98 99 100 101 102 103 c Tracers : 104 c ~~~~~~~ 100 ! variables added for CO2 condensation 101 ! ------------------------------------ 102 REAL hh !, zhcond(ngridmx,nlayermx) 103 ! REAL latcond,tcond1mb 104 ! REAL acond,bcond 105 ! SAVE acond,bcond 106 ! DATA latcond,tcond1mb/5.9e5,136.27/ 107 108 ! Tracers 109 ! ------- 105 110 INTEGER iq 106 111 REAL zq(ngridmx,nlayermx,nqmx) 107 112 REAL zq1temp(ngridmx) 108 REAL rho(ngridmx) ! nearsurface air density113 REAL rho(ngridmx) ! near-surface air density 109 114 REAL qsat(ngridmx) 110 115 DATA firstcall/.true./ 111 116 REAL kmixmin 112 117 113 c ** un petit test de coherence 114 c -------------------------- 118 ! Variables added for implicit latent heat inclusion 119 ! -------------------------------------------------- 120 real latconst, dqsat(ngridmx), qsat_temp1, qsat_temp2 121 real z1_Tdry(ngridmx), z2_Tdry(ngridmx) 122 real z1_T(ngridmx), z2_T(ngridmx) 123 real zb_T(ngridmx) 124 real zc_T(ngridmx,nlayermx) 125 real zd_T(ngridmx,nlayermx) 126 real lat1(ngridmx), lat2(ngridmx) 127 real surfh2otot 128 logical surffluxdiag 129 integer isub ! sub-loop for precision 130 131 integer ivap, iice ! also make liq for clarity on surface... 132 save ivap, iice 133 134 real, parameter :: sigma=5.67e-8 135 real, parameter :: karman=0.4 136 real cd0, roughratio 137 138 logical forceWC 139 real masse, Wtot, Wdiff 140 141 real dqsdif_total(ngrid) 142 real zq0(ngrid) 143 144 forceWC=.true. 145 ! forceWC=.false. 146 147 148 ! Coherence test 149 ! -------------- 115 150 116 151 IF (firstcall) THEN 117 IF(ngrid.NE.ngridmx) THEN 118 PRINT*,'STOP dans vdifc' 119 PRINT*,'probleme de dimensions :' 120 PRINT*,'ngrid =',ngrid 121 PRINT*,'ngridmx =',ngridmx 122 STOP 123 ENDIF 124 c To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal) 125 bcond=1./tcond1mb 126 acond=r/latcond 127 PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond 128 PRINT*,' acond,bcond',acond,bcond 129 130 firstcall=.false. 152 ! IF(ngrid.NE.ngridmx) THEN 153 ! PRINT*,'STOP dans vdifc' 154 ! PRINT*,'probleme de dimensions :' 155 ! PRINT*,'ngrid =',ngrid 156 ! PRINT*,'ngridmx =',ngridmx 157 ! STOP 158 ! ENDIF 159 ! To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal) 160 ! bcond=1./tcond1mb 161 ! acond=r/latcond 162 ! PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond 163 ! PRINT*,' acond,bcond',acond,bcond 164 165 if(water)then 166 ! iliq=igcm_h2o_vap 167 ivap=igcm_h2o_vap 168 iice=igcm_h2o_ice ! simply to make the code legible 169 ! to be generalised later 170 endif 171 172 if(ngridmx.ne.1)then 173 if(nonideal)then 174 print*,'Nonideal doesnt work yet in 3D!!!' 175 call abort 176 endif 177 endif 178 179 firstcall=.false. 131 180 ENDIF 132 181 133 134 135 136 137 c----------------------------------------------------------------------- 138 c 1. initialisation 139 c ----------------- 182 !----------------------------------------------------------------------- 183 ! 1. Initialisation 184 ! ----------------- 140 185 141 186 nlev=nlay+1 142 187 143 c ** calcul de rho*dz etdt*rho/dz=dt*rho**2 g/dp144 c avecrho=p/RT=p/ (R Theta) (p/ps)**kappa145 c----------------------------------------188 ! Calculate rho*dz and dt*rho/dz=dt*rho**2 g/dp 189 ! with rho=p/RT=p/ (R Theta) (p/ps)**kappa 190 ! --------------------------------------------- 146 191 147 192 DO ilay=1,nlay … … 151 196 ENDDO 152 197 153 zcst1=4.*g*ptimestep/( r*r)198 zcst1=4.*g*ptimestep/(R*R) 154 199 DO ilev=2,nlev-1 155 200 DO ig=1,ngrid 156 201 zb0(ig,ilev)=pplev(ig,ilev)* 157 s (pplev(ig,1)/pplev(ig,ilev))**rcp /158 s (ph(ig,ilev-1)+ph(ig,ilev))202 s (pplev(ig,1)/pplev(ig,ilev))**rcp / 203 s (ph(ig,ilev-1)+ph(ig,ilev)) 159 204 zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/ 160 s (pplay(ig,ilev-1)-pplay(ig,ilev)) 161 ENDDO 162 ENDDO 163 DO ig=1,ngrid 164 zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig)) 165 ENDDO 166 167 c ** diagnostique pour l'initialisation 168 c ---------------------------------- 169 170 IF(lecrit) THEN 171 ig=ngrid/2+1 172 PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz' 173 DO ilay=1,nlay 174 WRITE(*,'(6f11.5)') 175 s .01*pplay(ig,ilay),.001*pzlay(ig,ilay), 176 s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay) 177 ENDDO 178 PRINT*,'Pression (mbar) ,altitude (km),zb' 179 DO ilev=1,nlay 180 WRITE(*,'(3f15.7)') 181 s .01*pplev(ig,ilev),.001*pzlev(ig,ilev), 182 s zb0(ig,ilev) 183 ENDDO 184 ENDIF 185 186 c Potential Condensation temperature: 187 c ----------------------------------- 188 189 c if (callcond) then 190 c DO ilev=1,nlay 191 c DO ig=1,ngrid 192 c zhcond(ig,ilev) = 193 c & (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev) 194 c END DO 195 c END DO 196 c else 197 call zerophys(ngrid*nlay,zhcond) 198 c end if 199 200 201 c----------------------------------------------------------------------- 202 c 2. ajout des tendances physiques 203 c ----------------------------- 205 s (pplay(ig,ilev-1)-pplay(ig,ilev)) 206 ENDDO 207 ENDDO 208 DO ig=1,ngrid 209 zb0(ig,1)=ptimestep*pplev(ig,1)/(R*ptsrf(ig)) 210 ENDDO 211 212 dqsdif_total(:)=0.0 213 214 !----------------------------------------------------------------------- 215 ! 2. Add the physical tendencies computed so far 216 ! ---------------------------------------------- 204 217 205 218 DO ilev=1,nlay … … 208 221 zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep 209 222 zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 210 zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))211 223 ENDDO 212 224 ENDDO 213 225 if(tracer) then 214 DO iq =1, nq 215 DO ilev=1,nlay 216 DO ig=1,ngrid 217 zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep 218 ENDDO 219 ENDDO 220 ENDDO 226 DO iq =1, nq 227 DO ilev=1,nlay 228 DO ig=1,ngrid 229 zq(ig,ilev,iq)=pq(ig,ilev,iq) + 230 & pdqfi(ig,ilev,iq)*ptimestep 231 ENDDO 232 ENDDO 233 ENDDO 221 234 end if 222 235 223 c----------------------------------------------------------------------- 224 c 3. schema de turbulence 225 c -------------------- 226 227 c ** source d'energie cinetique turbulente a la surface 228 c (condition aux limites du schema de diffusion turbulente 229 c dans la couche limite 230 c --------------------- 231 232 CALL vdif_cd( ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph 233 & ,zcdv_true,zcdh_true) 234 DO ig=1,ngrid 235 zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) 236 zcdv(ig)=zcdv_true(ig)*sqrt(zu2) 237 zcdh(ig)=zcdh_true(ig)*sqrt(zu2) 238 ENDDO 239 240 c ** schema de diffusion turbulente dans la couche limite 241 c ---------------------------------------------------- 242 243 CALL vdif_kc(ptimestep,g,pzlev,pzlay 244 & ,pu,pv,ph,zcdv_true 245 & ,pq2,zkv,zkh) 246 247 248 c Adding eddy mixing to mimic 3D general circulation in 1D 249 c RW FF 2010 236 !----------------------------------------------------------------------- 237 ! 3. Turbulence scheme 238 ! -------------------- 239 ! 240 ! Source of turbulent kinetic energy at the surface 241 ! ------------------------------------------------- 242 ! Formula is Cd_0 = (karman / log[1+z1/z0])^2 243 244 DO ig=1,ngrid 245 roughratio = 1.E+0 + pzlay(ig,1)/pz0 246 cd0 = karman/log(roughratio) 247 cd0 = cd0*cd0 248 zcdv_true(ig) = cd0 249 zcdh_true(ig) = cd0 250 ENDDO 251 252 DO ig=1,ngrid 253 zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) 254 zcdv(ig)=zcdv_true(ig)*sqrt(zu2) 255 zcdh(ig)=zcdh_true(ig)*sqrt(zu2) 256 ENDDO 257 258 ! Turbulent diffusion coefficients in the boundary layer 259 ! ------------------------------------------------------ 260 261 call vdif_kc(ptimestep,g,pzlev,pzlay 262 & ,pu,pv,ph,zcdv_true 263 & ,pq2,zkv,zkh) 264 265 ! Adding eddy mixing to mimic 3D general circulation in 1D 266 ! R. Wordsworth & F. Forget (2010) 250 267 if ((ngrid.eq.1)) then 251 kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D 252 do ilev=1,nlay 253 do ig=1,ngrid 254 zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev)) 255 zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev)) 256 !zkh(ig,ilev) = kmixmin 257 !zkv(ig,ilev) = kmixmin 258 end do 259 end do 268 kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D 269 do ilev=1,nlay 270 do ig=1,ngrid 271 !zkh(ig,ilev) = 1.0 272 zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev)) 273 zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev)) 274 end do 275 end do 260 276 end if 261 277 262 263 c ** diagnostique pour le schema de turbulence 264 c ----------------------------------------- 265 266 IF(lecrit) THEN 267 PRINT* 268 PRINT*,'Diagnostic for the vertical turbulent mixing' 269 PRINT*,'Cd for momentum and potential temperature' 270 271 PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1) 272 PRINT*,'Mixing coefficient for momentum and pot.temp.' 273 DO ilev=1,nlay 274 PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev) 275 ENDDO 276 ENDIF 277 278 c----------------------------------------------------------------------- 279 c 4. inversion pour l'implicite sur u 280 c -------------------------------- 281 282 c ** l'equation est 283 c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1) 284 c avec 285 c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.) 286 c et 287 c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1) 288 c donc les entrees sont /zcdv/ pour la condition a la limite sol 289 c et /zkv/ = Ku 290 278 !----------------------------------------------------------------------- 279 ! 4. Implicit inversion of u 280 ! -------------------------- 281 282 ! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1) 283 ! avec 284 ! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.) 285 ! et 286 ! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1) 287 ! donc les entrees sont /zcdv/ pour la condition a la limite sol 288 ! et /zkv/ = Ku 289 291 290 CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2)) 292 291 CALL multipl(ngrid,zcdv,zb0,zb) … … 301 300 DO ig=1,ngrid 302 301 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 303 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))302 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 304 303 zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+ 305 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)304 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 306 305 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 307 306 ENDDO … … 317 316 ENDDO 318 317 319 c----------------------------------------------------------------------- 320 c 5. inversion pour l'implicite sur v 321 c -------------------------------- 322 323 c ** l'equation est 324 c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1) 325 c avec 326 c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.) 327 c et 328 c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1) 329 c donc les entrees sont /zcdv/ pour la condition a la limite sol 330 c et /zkv/ = Kv 318 !----------------------------------------------------------------------- 319 ! 5. Implicit inversion of v 320 ! -------------------------- 321 322 ! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1) 323 ! avec 324 ! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.) 325 ! et 326 ! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1) 327 ! donc les entrees sont /zcdv/ pour la condition a la limite sol 328 ! et /zkv/ = Kv 331 329 332 330 DO ig=1,ngrid … … 339 337 DO ig=1,ngrid 340 338 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 341 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))339 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 342 340 zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+ 343 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)341 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 344 342 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 345 343 ENDDO … … 355 353 ENDDO 356 354 357 c----------------------------------------------------------------------- 358 c 6. inversion pour l'implicite sur h sans oublier le couplage 359 c avec le sol (conduction) 360 c ------------------------ 361 362 c ** l'equation est 363 c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1) 364 c avec 365 c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.) 366 c et 367 c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1) 368 c donc les entrees sont /zcdh/ pour la condition de raccord au sol 369 c et /zkh/ = Kh 370 c ------------- 355 !---------------------------------------------------------------------------- 356 ! 6. Implicit inversion of h, not forgetting the coupling with the ground 357 358 ! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1) 359 ! avec 360 ! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.) 361 ! et 362 ! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1) 363 ! donc les entrees sont /zcdh/ pour la condition de raccord au sol 364 ! et /zkh/ = Kh 365 366 ! Using the wind modified by friction for lifting and sublimation 367 ! --------------------------------------------------------------- 368 DO ig=1,ngrid 369 zu2 = zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) 370 zcdv(ig) = zcdv_true(ig)*sqrt(zu2) 371 zcdh(ig) = zcdh_true(ig)*sqrt(zu2) 372 ENDDO 371 373 372 374 CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2)) … … 379 381 ENDDO 380 382 381 DO ilay=nlay-1, 1,-1383 DO ilay=nlay-1,2,-1 382 384 DO ig=1,ngrid 383 385 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 384 $zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))386 & zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 385 387 zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ 386 $zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)388 & zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 387 389 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 388 390 ENDDO 389 391 ENDDO 390 392 391 c ** calcul de (d Planck / dT) a la temperature d'interface 392 c ------------------------------------------------------ 393 394 z4st=4.*5.67e-8*ptimestep 393 DO ig=1,ngrid 394 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 395 & zb(ig,2)*(1.-zd(ig,2))) 396 zc(ig,1)=(za(ig,1)*zh(ig,1)+ 397 & zb(ig,2)*zc(ig,2))*z1(ig) 398 zd(ig,1)=zb(ig,1)*z1(ig) 399 ENDDO 400 401 ! Calculate (d Planck / dT) at the interface temperature 402 ! ------------------------------------------------------ 403 404 z4st=4.0*sigma*ptimestep 395 405 DO ig=1,ngrid 396 406 zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) 397 407 ENDDO 398 408 399 c ** calcul de la temperature_d'interface et de sa tendance. 400 c on ecrit que la somme des flux est nulle a l'interface 401 c a t + \delta t, 402 c c'est a dire le flux radiatif a {t + \delta t} 403 c + le flux turbulent a {t + \delta t} 404 c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1 405 c (notation K dt = /cpp*b/) 406 c + le flux dans le sol a t 407 c + l'evolution du flux dans le sol lorsque la temperature d'interface 408 c passe de sa valeur a t a sa valeur a {t + \delta t}. 409 c ---------------------------------------------------- 410 411 DO ig=1,ngrid 412 z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) 413 s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep 414 z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig) 415 ztsrf2(ig)=z1(ig)/z2(ig) 416 pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep 417 418 c Modif speciale CO2 condensation: 419 c tconds = 1./(bcond-acond*log(.0095*pplev(ig,1))) 420 c if ((callcond).and. 421 c & ((co2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then 422 c zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds 423 c else 424 zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) 425 c end if 426 ENDDO 427 428 c ** et a partir de la temperature au sol on remonte 429 c ----------------------------------------------- 430 431 DO ilay=2,nlay 432 DO ig=1,ngrid 433 hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif co2cond 434 zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh 435 ENDDO 436 ENDDO 437 438 439 c----------------------------------------------------------------------- 440 c TRACERS 441 c ------- 409 ! Calculate temperature tendency at the interface (dry case) 410 ! ---------------------------------------------------------- 411 ! Sum of fluxes at interface at time t + \delta t gives change in T: 412 ! radiative fluxes 413 ! turbulent convective (sensible) heat flux 414 ! flux (if any) from subsurface 415 416 if(.not.water) then 417 418 DO ig=1,ngrid 419 420 z1(ig) = pcapcal(ig)*ptsrf(ig) + cpp*zb(ig,1)*zc(ig,1) 421 & + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep 422 z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1)) 423 & +zdplanck(ig) 424 ztsrf2(ig) = z1(ig) / z2(ig) 425 pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep 426 zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig) 427 ENDDO 428 429 ! Recalculate temperature to top of atmosphere, starting from ground 430 ! ------------------------------------------------------------------ 431 432 DO ilay=2,nlay 433 DO ig=1,ngrid 434 hh = zh(ig,ilay-1) 435 zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh 436 ENDDO 437 ENDDO 438 439 endif ! not water 440 441 !----------------------------------------------------------------------- 442 ! TRACERS (no vapour) 443 ! ------- 442 444 443 445 if(tracer) then 444 446 445 c Using the wind modified by friction for lifting and sublimation 446 c ---------------------------------------------------------------- 447 DO ig=1,ngrid 448 zu2=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) 449 zcdv(ig)=zcdv_true(ig)*sqrt(zu2) 450 zcdh(ig)=zcdh_true(ig)*sqrt(zu2) 451 ENDDO 452 453 c Calcul du flux vertical au bas de la premiere couche (dust) : 454 c ----------------------------------------------------------- 455 do ig=1,ngridmx 456 rho(ig) = zb0(ig,1) /ptimestep 457 ! zb(ig,1) = 0. 458 end do 459 460 call zerophys(ngrid*nq,pdqsdif) 461 462 c OU calcul de la valeur de q a la surface (water) : 463 c ---------------------------------------- 464 if (water) then 465 !call watersat(ngridmx,ptsrf,pplev(1,1),qsat) 466 do ig=1,ngridmx 467 call watersat_2(ptsrf(ig),pplev(ig,1),qsat(ig)) 468 end do 469 end if 470 471 c Inversion pour l'implicite sur q 472 c -------------------------------- 473 do iq=1,nq 474 CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2)) 475 476 if ((water).and.(iq.eq.igcm_h2o_vap)) then 477 c This line is required to account for turbulent transport 478 c from surface (e.g. ice) to mid-layer of atmosphere: 479 CALL multipl(ngrid,zcdv,zb0,zb(1,1)) 480 CALL multipl(ngrid,dryness,zb(1,1),zb(1,1)) 481 else ! (re)-initialize zb(:,1) 482 zb(1:ngrid,1)=0 483 end if 484 485 DO ig=1,ngrid 486 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 487 zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) 488 zd(ig,nlay)=zb(ig,nlay)*z1(ig) 489 ENDDO 490 491 DO ilay=nlay-1,2,-1 447 ! Calculate vertical flux from the bottom to the first layer (dust) 448 ! ----------------------------------------------------------------- 449 do ig=1,ngridmx 450 rho(ig) = zb0(ig,1) /ptimestep 451 end do 452 453 call zerophys(ngrid*nq,pdqsdif) 454 455 ! Implicit inversion of q 456 ! ----------------------- 457 do iq=1,nq 458 459 if (iq.ne.igcm_h2o_vap) then 460 492 461 DO ig=1,ngrid 493 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 494 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 495 zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ 496 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 497 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 462 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 463 zcq(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) 464 zdq(ig,nlay)=zb(ig,nlay)*z1(ig) 465 ENDDO 466 467 DO ilay=nlay-1,2,-1 468 DO ig=1,ngrid 469 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 470 & zb(ig,ilay+1)*(1.-zdq(ig,ilay+1))) 471 zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ 472 & zb(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig) 473 zdq(ig,ilay)=zb(ig,ilay)*z1(ig) 474 ENDDO 498 475 ENDDO 499 ENDDO 500 501 if (water.and.(iq.eq.igcm_h2o_ice)) then 502 ! special case for water ice tracer: do not include 503 ! h2o ice tracer from surface (which is set when handling 504 ! h2o vapour case (see further down). 505 DO ig=1,ngrid 506 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 507 $ zb(ig,2)*(1.-zd(ig,2))) 508 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 509 $ zb(ig,2)*zc(ig,2)) *z1(ig) 510 ENDDO 511 else ! general case 512 DO ig=1,ngrid 513 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 514 $ zb(ig,2)*(1.-zd(ig,2))) 515 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 516 $ zb(ig,2)*zc(ig,2) + 517 $ (-pdqsdif(ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface 518 ENDDO 519 endif ! of if (water.and.(iq.eq.igcm_h2o_ice)) 520 521 IF ((water).and.(iq.eq.igcm_h2o_vap)) then 522 c Calculation for turbulent exchange with the surface (for ice) 523 DO ig=1,ngrid 524 zd(ig,1)=zb(ig,1)*z1(ig) 525 zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) 526 527 pdqsdif(ig,igcm_h2o_ice)=rho(ig)*dryness(ig)*zcdv(ig) 528 & *(zq1temp(ig)-qsat(ig)) 529 c write(*,*)'flux vers le sol=',pdqsdif(ig,nq) 530 END DO 531 532 DO ig=1,ngrid 533 if(.not.watercaptag(ig)) then 534 if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep) 535 & .gt.pqsurf(ig,igcm_h2o_ice)) then 536 c write(*,*)'on sublime plus que qsurf!' 537 pdqsdif(ig,igcm_h2o_ice)= 538 & -pqsurf(ig,igcm_h2o_ice)/ptimestep 539 c write(*,*)'flux vers le sol=',pdqsdif(ig,nq) 540 z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2))) 541 zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+ 542 $ zb(ig,2)*zc(ig,2) + 543 $ (-pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig) 544 zq1temp(ig)=zc(ig,1) 545 endif 546 endif ! if (.not.watercaptag(ig)) 547 c Starting upward calculations for water : 548 zq(ig,1,igcm_h2o_vap)=zq1temp(ig) 549 ENDDO ! of DO ig=1,ngrid 550 551 552 ! ! ADDED BY RW: Water latent heat effect 553 ! this doesnt work; causes instability. What do they do on earth? 554 ! DO ig=1,ngrid 555 ! pdtsrf(ig)= 556 ! & pdtsrf(ig)+RLVTT*pdqsdif(ig,igcm_h2o_ice)/pcapcal(ig) 557 ! ENDDO 558 ! print*,'Surface latent heat release in vdifc' 559 ! print*,'due to H2O NOT taken into account.' 560 561 ELSE 562 c Starting upward calculations for simple mixing of tracer (dust) 563 DO ig=1,ngrid 564 zq(ig,1,iq)=zc(ig,1) 565 ENDDO 566 END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap)) 567 568 DO ilay=2,nlay 569 DO ig=1,ngrid 570 zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) 571 ENDDO 572 ENDDO 573 enddo ! of do iq=1,nq 574 end if ! of if(tracer) 575 576 c----------------------------------------------------------------------- 577 c 8. calcul final des tendances de la diffusion verticale 578 c ---------------------------------------------------- 579 580 DO ilev = 1, nlay 581 DO ig=1,ngrid 582 pdudif(ig,ilev)=( zu(ig,ilev)- 583 $ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep 584 pdvdif(ig,ilev)=( zv(ig,ilev)- 585 $ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep 586 hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep , 587 $ zhcond(ig,ilev)) ! modif co2cond 476 477 if ((water).and.(iq.eq.iice)) then 478 ! special case for water ice tracer: do not include 479 ! h2o ice tracer from surface (which is set when handling 480 ! h2o vapour case (see further down). 481 ! zb(ig,1)=0 if iq ne ivap 482 DO ig=1,ngrid 483 z1(ig)=1./(za(ig,1)+ 484 & zb(ig,2)*(1.-zdq(ig,2))) 485 zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 486 & zb(ig,2)*zcq(ig,2))*z1(ig) 487 ENDDO 488 else ! general case 489 DO ig=1,ngrid 490 z1(ig)=1./(za(ig,1)+ 491 & zb(ig,2)*(1.-zdq(ig,2))) 492 zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 493 & zb(ig,2)*zcq(ig,2) 494 & +(-pdqsdif(ig,iq))*ptimestep)*z1(ig) 495 ! tracer flux from surface 496 ! currently pdqsdif always zero here, 497 ! so last line is superfluous 498 enddo 499 endif ! of if (water.and.(iq.eq.igcm_h2o_ice)) 500 501 502 ! Starting upward calculations for simple tracer mixing (e.g., dust) 503 do ig=1,ngrid 504 zq(ig,1,iq)=zcq(ig,1) 505 end do 506 507 do ilay=2,nlay 508 do ig=1,ngrid 509 zq(ig,ilay,iq)=zcq(ig,ilay)+ 510 $ zdq(ig,ilay)*zq(ig,ilay-1,iq) 511 end do 512 end do 513 514 endif ! if (iq.ne.igcm_h2o_vap) 515 516 ! Calculate temperature tendency including latent heat term 517 ! and assuming an infinite source of water on the ground 518 ! ------------------------------------------------------------------ 519 520 if (water.and.(iq.eq.igcm_h2o_vap)) then 521 522 ! compute evaporation efficiency 523 do ig = 1, ngrid 524 if(rnat(ig).eq.1)then 525 dryness(ig)=pqsurf(ig,ivap)+pqsurf(ig,iice) 526 dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol) 527 dryness(ig)=MAX(0.,dryness(ig)) 528 endif 529 enddo 530 531 do ig=1,ngrid 532 533 ! Calculate the value of qsat at the surface (water) 534 call watersat(ptsrf(ig),pplev(ig,1),qsat(ig)) 535 call watersat(ptsrf(ig)-0.0001,pplev(ig,1),qsat_temp1) 536 call watersat(ptsrf(ig)+0.0001,pplev(ig,1),qsat_temp2) 537 dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002 538 ! calculate dQsat / dT by finite differences 539 ! we cannot use the updated temperature value yet... 540 541 enddo 542 543 ! coefficients for q 544 545 do ig=1,ngrid 546 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 547 zcq(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) 548 zdq(ig,nlay)=zb(ig,nlay)*z1(ig) 549 enddo 550 551 do ilay=nlay-1,2,-1 552 do ig=1,ngrid 553 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 554 $ zb(ig,ilay+1)*(1.-zdq(ig,ilay+1))) 555 zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ 556 $ zb(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig) 557 zdq(ig,ilay)=zb(ig,ilay)*z1(ig) 558 enddo 559 enddo 560 561 do ig=1,ngrid 562 z1(ig)=1./(za(ig,1)+zb(ig,1)*dryness(ig)+ 563 $ zb(ig,2)*(1.-zdq(ig,2))) 564 zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 565 $ zb(ig,2)*zcq(ig,2))*z1(ig) 566 zdq(ig,1)=dryness(ig)*zb(ig,1)*z1(ig) 567 enddo 568 569 ! calculation of h0 and h1 570 do ig=1,ngrid 571 zdq0(ig) = dqsat(ig) 572 zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig) 573 574 z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zb(ig,1)*zc(ig,1) 575 & + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep 576 & + zb(ig,1)*dryness(ig)*RLVTT* 577 & ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1)) 578 579 z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1)) 580 & +zdplanck(ig) 581 & +zb(ig,1)*dryness(ig)*RLVTT*zdq0(ig)* 582 & (1.0-zdq(ig,1)) 583 584 ztsrf2(ig) = z1(ig) / z2(ig) 585 pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep 586 zh(ig,1) = zc(ig,1) + zd(ig,1)*ztsrf2(ig) 587 enddo 588 589 ! calculation of qs and q1 590 do ig=1,ngrid 591 zq0(ig) = zcq0(ig)+zdq0(ig)*ztsrf2(ig) 592 zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig) 593 enddo 594 595 ! calculation of evaporation 596 do ig=1,ngrid 597 evap(ig)= zb(ig,1)*dryness(ig)*(zq(ig,1,ivap)-zq0(ig)) 598 dqsdif_total(ig)=evap(ig) 599 enddo 600 601 ! recalculate temperature and q(vap) to top of atmosphere, starting from ground 602 do ilay=2,nlay 603 do ig=1,ngrid 604 zq(ig,ilay,iq)=zcq(ig,ilay) 605 & +zdq(ig,ilay)*zq(ig,ilay-1,iq) 606 zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1) 607 end do 608 end do 609 610 do ig=1,ngrid 611 612 ! -------------------------------------------------------------------------- 613 ! On the ocean, if T > 0 C then the vapour tendency must replace the ice one 614 ! The surface vapour tracer is actually liquid. To make things difficult. 615 616 if (rnat(ig).eq.0) then ! unfrozen ocean 617 618 pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep 619 pdqsdif(ig,iice)=0.0 620 621 622 elseif (rnat(ig).eq.1) then ! (continent) 623 624 ! -------------------------------------------------------- 625 ! Now check if we've taken too much water from the surface 626 ! This can only occur on the continent 627 628 ! If water is evaporating / subliming, we take it from ice before liquid 629 ! -- is this valid?? 630 if(dqsdif_total(ig).lt.0)then 631 pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep 632 pdqsdif(ig,iice)=max(-pqsurf(ig,iice)/ptimestep 633 & ,pdqsdif(ig,iice)) 634 endif 635 ! sublimation only greater than qsurf(ice) 636 ! ---------------------------------------- 637 ! we just convert some liquid to vapour too 638 ! if latent heats are the same, no big deal 639 if (-dqsdif_total(ig).gt.pqsurf(ig,iice))then 640 pdqsdif(ig,iice) = -pqsurf(ig,iice)/ptimestep ! removes all the ice! 641 pdqsdif(ig,ivap) = dqsdif_total(ig)/ptimestep 642 & - pdqsdif(ig,iice) ! take the remainder from the liquid instead 643 pdqsdif(ig,ivap) = max(-pqsurf(ig,ivap)/ptimestep 644 & ,pdqsdif(ig,ivap)) 645 endif 646 647 endif ! if (rnat.ne.1) 648 649 ! If water vapour is condensing, we must decide whether it forms ice or liquid. 650 if(dqsdif_total(ig).gt.0)then ! a bug was here! 651 if(ztsrf2(ig).gt.To)then 652 pdqsdif(ig,iice)=0.0 653 pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep 654 else 655 pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep 656 pdqsdif(ig,ivap)=0.0 657 endif 658 endif 659 660 end do ! of DO ig=1,ngrid 661 endif ! if (water et iq=ivap) 662 end do ! of do iq=1,nq 663 endif ! traceur 664 665 666 !----------------------------------------------------------------------- 667 ! 8. Final calculation of the vertical diffusion tendencies 668 ! ----------------------------------------------------------------- 669 670 do ilev = 1, nlay 671 do ig=1,ngrid 672 pdudif(ig,ilev)=(zu(ig,ilev)- 673 & (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep 674 pdvdif(ig,ilev)=(zv(ig,ilev)- 675 & (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep 676 hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 677 588 678 pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep 589 ENDDO 590 ENDDO 591 592 593 if (tracer) then 594 DO iq = 1, nq 595 DO ilev = 1, nlay 596 DO ig=1,ngrid 597 pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)- 598 $ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep 599 ENDDO 600 ENDDO 601 ENDDO 602 end if 603 604 c ** diagnostique final 605 c ------------------ 606 607 IF(lecrit) THEN 608 PRINT*,'In vdif' 609 PRINT*,'Ts (t) and Ts (t+st)' 610 WRITE(*,'(a10,3a15)') 611 s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)' 612 PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1) 613 DO ilev=1,nlay 614 WRITE(*,'(4f15.7)') 615 s ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev), 616 s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev) 617 618 ENDDO 619 ENDIF 620 621 RETURN 622 END 679 enddo 680 enddo 681 682 if (tracer) then 683 do iq = 1, nq 684 do ilev = 1, nlay 685 do ig=1,ngrid 686 pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)- 687 & (pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/ 688 & ptimestep 689 enddo 690 enddo 691 enddo 692 693 if(water.and.forceWC)then ! force water conservation in model 694 ! we calculate the difference and add it to the ground 695 ! this is ugly and should be improved in the future 696 do ig=1,ngrid 697 Wtot=0.0 698 do ilay=1,nlay 699 masse = (pplev(ig,ilay) - pplev(ig,ilay+1))/g 700 ! Wtot=Wtot+masse*(zq(ig,ilay,iice)- 701 ! & (pq(ig,ilay,iice)+pdqfi(ig,ilay,iice)*ptimestep)) 702 Wtot=Wtot+masse*(zq(ig,ilay,ivap)- 703 & (pq(ig,ilay,ivap)+pdqfi(ig,ilay,ivap)*ptimestep)) 704 enddo 705 Wdiff=Wtot/ptimestep+pdqsdif(ig,ivap)+pdqsdif(ig,iice) 706 707 if(ztsrf2(ig).gt.To)then 708 pdqsdif(ig,ivap)=pdqsdif(ig,ivap)-Wdiff 709 else 710 pdqsdif(ig,iice)=pdqsdif(ig,iice)-Wdiff 711 endif 712 enddo 713 714 endif 715 716 endif 717 718 if(water)then 719 call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness) 720 endif 721 722 return 723 end -
trunk/LMDZ.GENERIC/libf/phystd/watercap.h
r135 r253 1 c-----------------------------------------------------------------------2 cINCLUDE 'watercap.h'1 !----------------------------------------------------------------------- 2 ! INCLUDE 'watercap.h' 3 3 4 4 logical watercaptag(ngridmx) -
trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90
r135 r253 6 6 real, parameter :: To = 273.16 7 7 real, parameter :: mH2O = 18.01528 8 real, parameter :: RLVTT = 2.5008E+6 ! Latent heat of sublimation (J kg-1) 8 9 ! benjamin additions 10 real, parameter :: RLVTT = 2.257E+6 ! Latent heat of sublimation (J kg-1) 11 real, parameter :: RLSTT = 2.257E+6 ! 2.591E+6 in reality ! Latent heat of sublimation (J kg-1) 12 !real, parameter :: RLVTT = 0.0 13 !real, parameter :: RLSTT = 0.0 14 15 real, parameter :: RLFTT = 3.334E+5 ! Latent heat of fusion (J kg-1) 16 !real, parameter :: RLFTT = 0.0 17 real, parameter :: rhowater = 1.0E+3 ! mass of water (kg/m^3) 18 real, parameter :: mx_eau_sol = 150 ! mass of water (kg/m^2) 9 19 10 20 ! This was the old Martian latent heat version: -
trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F
r135 r253 6 6 ! 0d (pour un scalaire qui ne depend que du temps : ex : la longitude 7 7 ! solaire) 8 ! (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)8 ! (ou encore 1d, dans le cas de rcm1d, pour sortir une colonne) 9 9 ! La periode d'ecriture est donnee par 10 10 ! "ecritphy " regle dans le fichier de controle de run : run.def … … 149 149 150 150 if (ngridmx.eq.1) then 151 ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme151 ! in rcm1d, for the 1d version of the GCM, iphysiq and irythme 152 152 ! are undefined; so set them to 1 153 153 iphysiq=1
Note: See TracChangeset
for help on using the changeset viewer.