[1969] | 1 | MODULE vdifc_mod |
---|
| 2 | |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | |
---|
| 5 | CONTAINS |
---|
| 6 | |
---|
[38] | 7 | SUBROUTINE vdifc(ngrid,nlay,nq,co2ice,ppopsk, |
---|
| 8 | $ ptimestep,pcapcal,lecrit, |
---|
| 9 | $ pplay,pplev,pzlay,pzlev,pz0, |
---|
| 10 | $ pu,pv,ph,pq,ptsrf,pemis,pqsurf, |
---|
| 11 | $ pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf, |
---|
| 12 | $ pdudif,pdvdif,pdhdif,pdtsrf,pq2, |
---|
[660] | 13 | $ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true, |
---|
[1996] | 14 | $ hfmax,pcondicea_co2microp,sensibFlux, |
---|
[2260] | 15 | $ dustliftday,local_time,watercap, dwatercap_dif) |
---|
[1236] | 16 | |
---|
[1036] | 17 | use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number, |
---|
| 18 | & igcm_dust_submicron, igcm_h2o_vap, |
---|
[1974] | 19 | & igcm_h2o_ice, alpha_lift, |
---|
[2312] | 20 | & igcm_hdo_vap, igcm_hdo_ice, |
---|
[1974] | 21 | & igcm_stormdust_mass, igcm_stormdust_number |
---|
[1224] | 22 | use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness |
---|
[1969] | 23 | USE comcstfi_h, ONLY: cpp, r, rcp, g |
---|
[1996] | 24 | use watersat_mod, only: watersat |
---|
[1242] | 25 | use turb_mod, only: turb_resolved, ustar, tstar |
---|
[2160] | 26 | use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol |
---|
[2312] | 27 | use hdo_surfex_mod, only: hdo_surfex |
---|
[2515] | 28 | c use geometry_mod, only: longitude_deg,latitude_deg ! Joseph |
---|
[2409] | 29 | use dust_param_mod, only: doubleq, submicron, lifting |
---|
[2312] | 30 | |
---|
[38] | 31 | IMPLICIT NONE |
---|
| 32 | |
---|
| 33 | c======================================================================= |
---|
| 34 | c |
---|
| 35 | c subject: |
---|
| 36 | c -------- |
---|
| 37 | c Turbulent diffusion (mixing) for potential T, U, V and tracer |
---|
| 38 | c |
---|
| 39 | c Shema implicite |
---|
| 40 | c On commence par rajouter au variables x la tendance physique |
---|
| 41 | c et on resoult en fait: |
---|
| 42 | c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1) |
---|
| 43 | c |
---|
| 44 | c author: |
---|
| 45 | c ------ |
---|
| 46 | c Hourdin/Forget/Fournier |
---|
| 47 | c======================================================================= |
---|
| 48 | |
---|
| 49 | c----------------------------------------------------------------------- |
---|
| 50 | c declarations: |
---|
| 51 | c ------------- |
---|
| 52 | |
---|
[1944] | 53 | include "callkeys.h" |
---|
| 54 | include "microphys.h" |
---|
[38] | 55 | |
---|
| 56 | c |
---|
| 57 | c arguments: |
---|
| 58 | c ---------- |
---|
| 59 | |
---|
[1036] | 60 | INTEGER,INTENT(IN) :: ngrid,nlay |
---|
| 61 | REAL,INTENT(IN) :: ptimestep |
---|
| 62 | REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1) |
---|
| 63 | REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1) |
---|
| 64 | REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) |
---|
| 65 | REAL,INTENT(IN) :: ph(ngrid,nlay) |
---|
| 66 | REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid) |
---|
| 67 | REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay) |
---|
| 68 | REAL,INTENT(IN) :: pdhfi(ngrid,nlay) |
---|
| 69 | REAL,INTENT(IN) :: pfluxsrf(ngrid) |
---|
| 70 | REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay) |
---|
| 71 | REAL,INTENT(OUT) :: pdtsrf(ngrid),pdhdif(ngrid,nlay) |
---|
[1130] | 72 | REAL,INTENT(IN) :: pcapcal(ngrid) |
---|
| 73 | REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1) |
---|
[38] | 74 | |
---|
| 75 | c Argument added for condensation: |
---|
[1036] | 76 | REAL,INTENT(IN) :: co2ice (ngrid), ppopsk(ngrid,nlay) |
---|
| 77 | logical,INTENT(IN) :: lecrit |
---|
[1996] | 78 | REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlay)! tendency due to CO2 condensation (kg/kg.s-1) |
---|
| 79 | |
---|
[1036] | 80 | REAL,INTENT(IN) :: pz0(ngrid) ! surface roughness length (m) |
---|
[224] | 81 | |
---|
[256] | 82 | c Argument added to account for subgrid gustiness : |
---|
| 83 | |
---|
[1944] | 84 | REAL,INTENT(IN) :: wstar(ngrid), hfmax(ngrid)!, zi(ngrid) |
---|
[256] | 85 | |
---|
[38] | 86 | c Traceurs : |
---|
[1036] | 87 | integer,intent(in) :: nq |
---|
| 88 | REAL,INTENT(IN) :: pqsurf(ngrid,nq) |
---|
[2515] | 89 | REAL :: zqsurf(ngrid) ! temporary water tracer |
---|
[1036] | 90 | real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) |
---|
| 91 | real,intent(out) :: pdqdif(ngrid,nlay,nq) |
---|
[1974] | 92 | real,intent(out) :: pdqsdif(ngrid,nq) |
---|
| 93 | REAL,INTENT(in) :: dustliftday(ngrid) |
---|
| 94 | REAL,INTENT(in) :: local_time(ngrid) |
---|
[38] | 95 | |
---|
| 96 | c local: |
---|
| 97 | c ------ |
---|
| 98 | |
---|
[1130] | 99 | REAL :: pt(ngrid,nlay) |
---|
[473] | 100 | |
---|
[38] | 101 | INTEGER ilev,ig,ilay,nlev |
---|
| 102 | |
---|
[1047] | 103 | REAL z4st,zdplanck(ngrid) |
---|
| 104 | REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1) |
---|
| 105 | REAL zkq(ngrid,nlay+1) |
---|
| 106 | REAL zcdv(ngrid),zcdh(ngrid) |
---|
| 107 | REAL zcdv_true(ngrid),zcdh_true(ngrid) ! drag coeff are used by the LES to recompute u* and hfx |
---|
| 108 | REAL zu(ngrid,nlay),zv(ngrid,nlay) |
---|
| 109 | REAL zh(ngrid,nlay) |
---|
| 110 | REAL ztsrf2(ngrid) |
---|
| 111 | REAL z1(ngrid),z2(ngrid) |
---|
| 112 | REAL za(ngrid,nlay),zb(ngrid,nlay) |
---|
| 113 | REAL zb0(ngrid,nlay) |
---|
| 114 | REAL zc(ngrid,nlay),zd(ngrid,nlay) |
---|
[38] | 115 | REAL zcst1 |
---|
[1047] | 116 | REAL zu2(ngrid) |
---|
[38] | 117 | |
---|
| 118 | EXTERNAL SSUM,SCOPY |
---|
| 119 | REAL SSUM |
---|
[1036] | 120 | LOGICAL,SAVE :: firstcall=.true. |
---|
[38] | 121 | |
---|
[626] | 122 | |
---|
[38] | 123 | c variable added for CO2 condensation: |
---|
| 124 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
[1047] | 125 | REAL hh , zhcond(ngrid,nlay) |
---|
| 126 | REAL,PARAMETER :: latcond=5.9e5 |
---|
| 127 | REAL,PARAMETER :: tcond1mb=136.27 |
---|
| 128 | REAL,SAVE :: acond,bcond |
---|
[669] | 129 | |
---|
[2515] | 130 | c Subtimestep & implicit treatment of water vapor |
---|
| 131 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 132 | REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice |
---|
| 133 | REAL ztsrf(ngrid) ! temporary surface temperature in tsub |
---|
| 134 | REAL zdtsrf(ngrid) ! surface temperature tendancy in tsub |
---|
| 135 | |
---|
[2179] | 136 | c For latent heat release from ground water ice sublimation |
---|
[2515] | 137 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 138 | REAL tsrf_lh(ngrid) ! temporary surface temperature with lh effect |
---|
[2179] | 139 | REAL lh ! latent heat, formulation given in the Technical Document: |
---|
| 140 | ! "Modeling water ice sublimation under Phoenix-like conditions", Montmessin et al. 2004 |
---|
[38] | 141 | |
---|
| 142 | c Tracers : |
---|
| 143 | c ~~~~~~~ |
---|
| 144 | INTEGER iq |
---|
[1047] | 145 | REAL zq(ngrid,nlay,nq) |
---|
| 146 | REAL zq1temp(ngrid) |
---|
| 147 | REAL rho(ngrid) ! near surface air density |
---|
| 148 | REAL qsat(ngrid) |
---|
[38] | 149 | |
---|
[2312] | 150 | REAL hdoflux(ngrid) ! value of vapour flux of HDO |
---|
| 151 | REAL h2oflux(ngrid) ! value of vapour flux of H2O |
---|
| 152 | REAL old_h2o_vap(ngrid) ! traceur d'eau avant traitement |
---|
| 153 | |
---|
[38] | 154 | REAL kmixmin |
---|
| 155 | |
---|
[2260] | 156 | c Argument added for surface water ice budget: |
---|
| 157 | REAL,INTENT(IN) :: watercap(ngrid) |
---|
| 158 | REAL,INTENT(OUT) :: dwatercap_dif(ngrid) |
---|
| 159 | |
---|
[2515] | 160 | c Subtimestep to compute h2o latent heat flux: |
---|
| 161 | REAL :: dtmax = 0.5 ! subtimestep temp criterion |
---|
| 162 | INTEGER tsub ! adaptative subtimestep (seconds) |
---|
| 163 | REAL subtimestep !ptimestep/nsubtimestep |
---|
| 164 | INTEGER nsubtimestep(ngrid) ! number of subtimestep (int) |
---|
| 165 | |
---|
[473] | 166 | c Mass-variation scheme : |
---|
| 167 | c ~~~~~~~ |
---|
| 168 | |
---|
| 169 | INTEGER j,l |
---|
[1047] | 170 | REAL zcondicea(ngrid,nlay) |
---|
| 171 | REAL zt(ngrid,nlay),ztcond(ngrid,nlay+1) |
---|
| 172 | REAL betam(ngrid,nlay),dmice(ngrid,nlay) |
---|
| 173 | REAL pdtc(ngrid,nlay) |
---|
| 174 | REAL zhs(ngrid,nlay) |
---|
| 175 | REAL,SAVE :: ccond |
---|
[473] | 176 | |
---|
| 177 | c Theta_m formulation for mass-variation scheme : |
---|
| 178 | c ~~~~~~~ |
---|
| 179 | |
---|
[1047] | 180 | INTEGER,SAVE :: ico2 |
---|
| 181 | INTEGER llnt(ngrid) |
---|
| 182 | REAL,SAVE :: m_co2, m_noco2, A , B |
---|
| 183 | REAL vmr_co2(ngrid,nlay) |
---|
[473] | 184 | REAL qco2,mmean |
---|
| 185 | |
---|
[1047] | 186 | REAL,INTENT(OUT) :: sensibFlux(ngrid) |
---|
[660] | 187 | |
---|
[2312] | 188 | !!MARGAUX |
---|
| 189 | REAL DoH_vap(ngrid,nlay) |
---|
| 190 | |
---|
[38] | 191 | c ** un petit test de coherence |
---|
| 192 | c -------------------------- |
---|
| 193 | |
---|
[1779] | 194 | ! AS: OK firstcall absolute |
---|
[38] | 195 | IF (firstcall) THEN |
---|
| 196 | c To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal) |
---|
| 197 | bcond=1./tcond1mb |
---|
| 198 | acond=r/latcond |
---|
[473] | 199 | ccond=cpp/(g*latcond) |
---|
[38] | 200 | PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond |
---|
[473] | 201 | PRINT*,' acond,bcond,ccond',acond,bcond,ccond |
---|
[38] | 202 | |
---|
[473] | 203 | |
---|
| 204 | ico2=0 |
---|
| 205 | |
---|
| 206 | if (tracer) then |
---|
| 207 | c Prepare Special treatment if one of the tracer is CO2 gas |
---|
[1036] | 208 | do iq=1,nq |
---|
[473] | 209 | if (noms(iq).eq."co2") then |
---|
| 210 | ico2=iq |
---|
| 211 | m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) |
---|
| 212 | m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) |
---|
| 213 | c Compute A and B coefficient use to compute |
---|
| 214 | c mean molecular mass Mair defined by |
---|
| 215 | c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2 |
---|
| 216 | c 1/Mair = A*q(ico2) + B |
---|
| 217 | A =(1/m_co2 - 1/m_noco2) |
---|
| 218 | B=1/m_noco2 |
---|
| 219 | endif |
---|
| 220 | enddo |
---|
| 221 | end if |
---|
| 222 | |
---|
[38] | 223 | firstcall=.false. |
---|
| 224 | ENDIF |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | c----------------------------------------------------------------------- |
---|
| 228 | c 1. initialisation |
---|
| 229 | c ----------------- |
---|
| 230 | |
---|
| 231 | nlev=nlay+1 |
---|
| 232 | |
---|
[1035] | 233 | ! initialize output tendencies to zero: |
---|
| 234 | pdudif(1:ngrid,1:nlay)=0 |
---|
| 235 | pdvdif(1:ngrid,1:nlay)=0 |
---|
| 236 | pdhdif(1:ngrid,1:nlay)=0 |
---|
| 237 | pdtsrf(1:ngrid)=0 |
---|
[2515] | 238 | zdtsrf(1:ngrid)=0 |
---|
[1035] | 239 | pdqdif(1:ngrid,1:nlay,1:nq)=0 |
---|
| 240 | pdqsdif(1:ngrid,1:nq)=0 |
---|
[2515] | 241 | zdqsdif(1:ngrid)=0 |
---|
[2260] | 242 | dwatercap_dif(1:ngrid)=0 |
---|
[1035] | 243 | |
---|
[38] | 244 | c ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp |
---|
| 245 | c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa |
---|
| 246 | c ---------------------------------------- |
---|
| 247 | |
---|
| 248 | DO ilay=1,nlay |
---|
| 249 | DO ig=1,ngrid |
---|
| 250 | za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g |
---|
[473] | 251 | ! Mass variation scheme: |
---|
| 252 | betam(ig,ilay)=-za(ig,ilay)*latcond/(cpp*ppopsk(ig,ilay)) |
---|
[38] | 253 | ENDDO |
---|
| 254 | ENDDO |
---|
| 255 | |
---|
| 256 | zcst1=4.*g*ptimestep/(r*r) |
---|
| 257 | DO ilev=2,nlev-1 |
---|
| 258 | DO ig=1,ngrid |
---|
| 259 | zb0(ig,ilev)=pplev(ig,ilev)* |
---|
| 260 | s (pplev(ig,1)/pplev(ig,ilev))**rcp / |
---|
| 261 | s (ph(ig,ilev-1)+ph(ig,ilev)) |
---|
| 262 | zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/ |
---|
| 263 | s (pplay(ig,ilev-1)-pplay(ig,ilev)) |
---|
| 264 | ENDDO |
---|
| 265 | ENDDO |
---|
| 266 | DO ig=1,ngrid |
---|
[473] | 267 | zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig)) |
---|
[38] | 268 | ENDDO |
---|
| 269 | |
---|
| 270 | c ** diagnostique pour l'initialisation |
---|
| 271 | c ---------------------------------- |
---|
| 272 | |
---|
| 273 | IF(lecrit) THEN |
---|
| 274 | ig=ngrid/2+1 |
---|
| 275 | PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz' |
---|
| 276 | DO ilay=1,nlay |
---|
| 277 | WRITE(*,'(6f11.5)') |
---|
| 278 | s .01*pplay(ig,ilay),.001*pzlay(ig,ilay), |
---|
| 279 | s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay) |
---|
| 280 | ENDDO |
---|
| 281 | PRINT*,'Pression (mbar) ,altitude (km),zb' |
---|
| 282 | DO ilev=1,nlay |
---|
| 283 | WRITE(*,'(3f15.7)') |
---|
| 284 | s .01*pplev(ig,ilev),.001*pzlev(ig,ilev), |
---|
| 285 | s zb0(ig,ilev) |
---|
| 286 | ENDDO |
---|
| 287 | ENDIF |
---|
| 288 | |
---|
[473] | 289 | c ----------------------------------- |
---|
[38] | 290 | c Potential Condensation temperature: |
---|
| 291 | c ----------------------------------- |
---|
| 292 | |
---|
[473] | 293 | c Compute CO2 Volume mixing ratio |
---|
| 294 | c ------------------------------- |
---|
| 295 | if (callcond.and.(ico2.ne.0)) then |
---|
| 296 | DO ilev=1,nlay |
---|
| 297 | DO ig=1,ngrid |
---|
[529] | 298 | qco2=MAX(1.E-30 |
---|
| 299 | & ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep) |
---|
[473] | 300 | c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2) |
---|
| 301 | mmean=1/(A*qco2 +B) |
---|
| 302 | vmr_co2(ig,ilev) = qco2*mmean/m_co2 |
---|
| 303 | ENDDO |
---|
| 304 | ENDDO |
---|
| 305 | else |
---|
| 306 | DO ilev=1,nlay |
---|
| 307 | DO ig=1,ngrid |
---|
| 308 | vmr_co2(ig,ilev)=0.95 |
---|
| 309 | ENDDO |
---|
| 310 | ENDDO |
---|
| 311 | end if |
---|
[38] | 312 | |
---|
[473] | 313 | c forecast of atmospheric temperature zt and frost temperature ztcond |
---|
| 314 | c -------------------------------------------------------------------- |
---|
[38] | 315 | |
---|
[473] | 316 | if (callcond) then |
---|
| 317 | DO ilev=1,nlay |
---|
| 318 | DO ig=1,ngrid |
---|
[884] | 319 | ztcond(ig,ilev)= |
---|
| 320 | & 1./(bcond-acond*log(.01*vmr_co2(ig,ilev)*pplay(ig,ilev))) |
---|
| 321 | if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica |
---|
[473] | 322 | ! zhcond(ig,ilev) = |
---|
| 323 | ! & (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev) |
---|
| 324 | zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev) |
---|
| 325 | END DO |
---|
| 326 | END DO |
---|
[884] | 327 | ztcond(:,nlay+1)=ztcond(:,nlay) |
---|
[473] | 328 | else |
---|
[884] | 329 | zhcond(:,:) = 0 |
---|
| 330 | ztcond(:,:) = 0 |
---|
[473] | 331 | end if |
---|
| 332 | |
---|
| 333 | |
---|
[38] | 334 | c----------------------------------------------------------------------- |
---|
| 335 | c 2. ajout des tendances physiques |
---|
| 336 | c ----------------------------- |
---|
| 337 | |
---|
| 338 | DO ilev=1,nlay |
---|
| 339 | DO ig=1,ngrid |
---|
| 340 | zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep |
---|
| 341 | zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep |
---|
| 342 | zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep |
---|
[473] | 343 | ! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev)) |
---|
[38] | 344 | ENDDO |
---|
| 345 | ENDDO |
---|
| 346 | if(tracer) then |
---|
| 347 | DO iq =1, nq |
---|
| 348 | DO ilev=1,nlay |
---|
| 349 | DO ig=1,ngrid |
---|
| 350 | zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep |
---|
| 351 | ENDDO |
---|
| 352 | ENDDO |
---|
| 353 | ENDDO |
---|
| 354 | end if |
---|
| 355 | |
---|
| 356 | c----------------------------------------------------------------------- |
---|
| 357 | c 3. schema de turbulence |
---|
| 358 | c -------------------- |
---|
| 359 | |
---|
| 360 | c ** source d'energie cinetique turbulente a la surface |
---|
| 361 | c (condition aux limites du schema de diffusion turbulente |
---|
| 362 | c dans la couche limite |
---|
| 363 | c --------------------- |
---|
| 364 | |
---|
[499] | 365 | CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph |
---|
[1236] | 366 | & ,zcdv_true,zcdh_true) |
---|
[256] | 367 | |
---|
[290] | 368 | zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1) |
---|
[256] | 369 | |
---|
[291] | 370 | IF (callrichsl) THEN |
---|
[545] | 371 | zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+ |
---|
| 372 | & (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2) |
---|
| 373 | zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+ |
---|
| 374 | & (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2) |
---|
[496] | 375 | |
---|
[1242] | 376 | ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+ |
---|
[545] | 377 | & (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2) |
---|
| 378 | |
---|
[1242] | 379 | tstar(:)=0. |
---|
[545] | 380 | DO ig=1,ngrid |
---|
| 381 | IF (zcdh_true(ig) .ne. 0.) THEN ! When Cd=Ch=0, u*=t*=0 |
---|
[1242] | 382 | tstar(ig)=(ph(ig,1)-ptsrf(ig))*zcdh(ig)/ustar(ig) |
---|
[545] | 383 | ENDIF |
---|
| 384 | ENDDO |
---|
| 385 | |
---|
[284] | 386 | ELSE |
---|
[545] | 387 | zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance |
---|
| 388 | zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance |
---|
[1242] | 389 | ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)) |
---|
| 390 | tstar(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:)) |
---|
[290] | 391 | ENDIF |
---|
[38] | 392 | |
---|
[529] | 393 | ! Some usefull diagnostics for the new surface layer parametrization : |
---|
[290] | 394 | |
---|
[1130] | 395 | ! call WRITEDIAGFI(ngrid,'vdifc_zcdv_true', |
---|
[256] | 396 | ! & 'momentum drag','no units', |
---|
| 397 | ! & 2,zcdv_true) |
---|
[1130] | 398 | ! call WRITEDIAGFI(ngrid,'vdifc_zcdh_true', |
---|
[256] | 399 | ! & 'heat drag','no units', |
---|
| 400 | ! & 2,zcdh_true) |
---|
[1130] | 401 | ! call WRITEDIAGFI(ngrid,'vdifc_ust', |
---|
[473] | 402 | ! & 'friction velocity','m/s',2,ust) |
---|
[1130] | 403 | ! call WRITEDIAGFI(ngrid,'vdifc_tst', |
---|
[473] | 404 | ! & 'friction temperature','K',2,tst) |
---|
[1130] | 405 | ! call WRITEDIAGFI(ngrid,'vdifc_zcdv', |
---|
[268] | 406 | ! & 'aerodyn momentum conductance','m/s', |
---|
[256] | 407 | ! & 2,zcdv) |
---|
[1130] | 408 | ! call WRITEDIAGFI(ngrid,'vdifc_zcdh', |
---|
[268] | 409 | ! & 'aerodyn heat conductance','m/s', |
---|
[256] | 410 | ! & 2,zcdh) |
---|
| 411 | |
---|
[38] | 412 | c ** schema de diffusion turbulente dans la couche limite |
---|
| 413 | c ---------------------------------------------------- |
---|
[1240] | 414 | IF (.not. callyamada4) THEN |
---|
[529] | 415 | |
---|
[1130] | 416 | CALL vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay |
---|
[38] | 417 | & ,pu,pv,ph,zcdv_true |
---|
[325] | 418 | & ,pq2,zkv,zkh,zq) |
---|
[529] | 419 | |
---|
[544] | 420 | ELSE |
---|
[38] | 421 | |
---|
[555] | 422 | pt(:,:)=ph(:,:)*ppopsk(:,:) |
---|
[1036] | 423 | CALL yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt |
---|
[1242] | 424 | s ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ustar |
---|
[652] | 425 | s ,9) |
---|
[544] | 426 | ENDIF |
---|
| 427 | |
---|
[38] | 428 | if ((doubleq).and.(ngrid.eq.1)) then |
---|
| 429 | kmixmin = 80. !80.! minimum eddy mix coeff in 1D |
---|
| 430 | do ilev=1,nlay |
---|
| 431 | do ig=1,ngrid |
---|
| 432 | zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev)) |
---|
| 433 | zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev)) |
---|
| 434 | end do |
---|
| 435 | end do |
---|
| 436 | end if |
---|
| 437 | |
---|
| 438 | c ** diagnostique pour le schema de turbulence |
---|
| 439 | c ----------------------------------------- |
---|
| 440 | |
---|
| 441 | IF(lecrit) THEN |
---|
| 442 | PRINT* |
---|
| 443 | PRINT*,'Diagnostic for the vertical turbulent mixing' |
---|
| 444 | PRINT*,'Cd for momentum and potential temperature' |
---|
| 445 | |
---|
| 446 | PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1) |
---|
| 447 | PRINT*,'Mixing coefficient for momentum and pot.temp.' |
---|
| 448 | DO ilev=1,nlay |
---|
| 449 | PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev) |
---|
| 450 | ENDDO |
---|
| 451 | ENDIF |
---|
| 452 | |
---|
| 453 | |
---|
| 454 | |
---|
| 455 | |
---|
| 456 | c----------------------------------------------------------------------- |
---|
| 457 | c 4. inversion pour l'implicite sur u |
---|
| 458 | c -------------------------------- |
---|
| 459 | |
---|
| 460 | c ** l'equation est |
---|
| 461 | c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1) |
---|
| 462 | c avec |
---|
| 463 | c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.) |
---|
| 464 | c et |
---|
| 465 | c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1) |
---|
| 466 | c donc les entrees sont /zcdv/ pour la condition a la limite sol |
---|
| 467 | c et /zkv/ = Ku |
---|
[2312] | 468 | |
---|
[2274] | 469 | zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) |
---|
| 470 | zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1) |
---|
[38] | 471 | |
---|
| 472 | DO ig=1,ngrid |
---|
| 473 | z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) |
---|
| 474 | zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig) |
---|
| 475 | zd(ig,nlay)=zb(ig,nlay)*z1(ig) |
---|
| 476 | ENDDO |
---|
| 477 | |
---|
| 478 | DO ilay=nlay-1,1,-1 |
---|
| 479 | DO ig=1,ngrid |
---|
| 480 | z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ |
---|
| 481 | $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) |
---|
| 482 | zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+ |
---|
| 483 | $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) |
---|
| 484 | zd(ig,ilay)=zb(ig,ilay)*z1(ig) |
---|
| 485 | ENDDO |
---|
| 486 | ENDDO |
---|
| 487 | |
---|
| 488 | DO ig=1,ngrid |
---|
| 489 | zu(ig,1)=zc(ig,1) |
---|
| 490 | ENDDO |
---|
| 491 | DO ilay=2,nlay |
---|
| 492 | DO ig=1,ngrid |
---|
| 493 | zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1) |
---|
| 494 | ENDDO |
---|
| 495 | ENDDO |
---|
| 496 | |
---|
| 497 | |
---|
| 498 | |
---|
| 499 | |
---|
| 500 | |
---|
| 501 | c----------------------------------------------------------------------- |
---|
| 502 | c 5. inversion pour l'implicite sur v |
---|
| 503 | c -------------------------------- |
---|
| 504 | |
---|
| 505 | c ** l'equation est |
---|
| 506 | c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1) |
---|
| 507 | c avec |
---|
| 508 | c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.) |
---|
| 509 | c et |
---|
| 510 | c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1) |
---|
| 511 | c donc les entrees sont /zcdv/ pour la condition a la limite sol |
---|
| 512 | c et /zkv/ = Kv |
---|
| 513 | |
---|
| 514 | DO ig=1,ngrid |
---|
| 515 | z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) |
---|
| 516 | zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig) |
---|
| 517 | zd(ig,nlay)=zb(ig,nlay)*z1(ig) |
---|
| 518 | ENDDO |
---|
| 519 | |
---|
| 520 | DO ilay=nlay-1,1,-1 |
---|
| 521 | DO ig=1,ngrid |
---|
| 522 | z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ |
---|
| 523 | $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) |
---|
| 524 | zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+ |
---|
| 525 | $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) |
---|
| 526 | zd(ig,ilay)=zb(ig,ilay)*z1(ig) |
---|
| 527 | ENDDO |
---|
| 528 | ENDDO |
---|
| 529 | |
---|
| 530 | DO ig=1,ngrid |
---|
| 531 | zv(ig,1)=zc(ig,1) |
---|
| 532 | ENDDO |
---|
| 533 | DO ilay=2,nlay |
---|
| 534 | DO ig=1,ngrid |
---|
| 535 | zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1) |
---|
| 536 | ENDDO |
---|
| 537 | ENDDO |
---|
| 538 | |
---|
| 539 | |
---|
| 540 | |
---|
| 541 | |
---|
| 542 | |
---|
| 543 | c----------------------------------------------------------------------- |
---|
| 544 | c 6. inversion pour l'implicite sur h sans oublier le couplage |
---|
| 545 | c avec le sol (conduction) |
---|
| 546 | c ------------------------ |
---|
| 547 | |
---|
| 548 | c ** l'equation est |
---|
| 549 | c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1) |
---|
| 550 | c avec |
---|
| 551 | c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.) |
---|
| 552 | c et |
---|
| 553 | c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1) |
---|
| 554 | c donc les entrees sont /zcdh/ pour la condition de raccord au sol |
---|
| 555 | c et /zkh/ = Kh |
---|
| 556 | c ------------- |
---|
| 557 | |
---|
[473] | 558 | c Mass variation scheme: |
---|
[2274] | 559 | zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) |
---|
| 560 | zb(1:ngrid,1)=zcdh(1:ngrid)*zb0(1:ngrid,1) |
---|
[38] | 561 | |
---|
[473] | 562 | c on initialise dm c |
---|
| 563 | |
---|
| 564 | pdtc(:,:)=0. |
---|
| 565 | zt(:,:)=0. |
---|
| 566 | dmice(:,:)=0. |
---|
[38] | 567 | |
---|
| 568 | c ** calcul de (d Planck / dT) a la temperature d'interface |
---|
| 569 | c ------------------------------------------------------ |
---|
| 570 | |
---|
| 571 | z4st=4.*5.67e-8*ptimestep |
---|
[544] | 572 | IF (tke_heat_flux .eq. 0.) THEN |
---|
[38] | 573 | DO ig=1,ngrid |
---|
| 574 | zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) |
---|
| 575 | ENDDO |
---|
[544] | 576 | ELSE |
---|
| 577 | zdplanck(:)=0. |
---|
| 578 | ENDIF |
---|
[38] | 579 | |
---|
[473] | 580 | ! calcul de zc et zd pour la couche top en prenant en compte le terme |
---|
[2080] | 581 | ! de variation de masse (on fait une boucle pour que \E7a converge) |
---|
[473] | 582 | |
---|
| 583 | ! Identification des points de grilles qui ont besoin de la correction |
---|
| 584 | |
---|
| 585 | llnt(:)=1 |
---|
[1236] | 586 | IF (.not.turb_resolved) THEN |
---|
[884] | 587 | IF (callcond) THEN |
---|
| 588 | DO ig=1,ngrid |
---|
[473] | 589 | DO l=1,nlay |
---|
| 590 | if(zh(ig,l) .lt. zhcond(ig,l)) then |
---|
| 591 | llnt(ig)=300 |
---|
| 592 | ! 200 and 100 do not go beyond month 9 with normal dissipation |
---|
| 593 | goto 5 |
---|
| 594 | endif |
---|
| 595 | ENDDO |
---|
[884] | 596 | 5 continue |
---|
| 597 | ENDDO |
---|
| 598 | ENDIF |
---|
[473] | 599 | |
---|
[529] | 600 | ENDIF |
---|
| 601 | |
---|
[473] | 602 | DO ig=1,ngrid |
---|
| 603 | |
---|
| 604 | ! Initialization of z1 and zd, which do not depend on dmice |
---|
| 605 | |
---|
| 606 | z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) |
---|
| 607 | zd(ig,nlay)=zb(ig,nlay)*z1(ig) |
---|
| 608 | |
---|
| 609 | DO ilay=nlay-1,1,-1 |
---|
| 610 | z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ |
---|
| 611 | $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) |
---|
| 612 | zd(ig,ilay)=zb(ig,ilay)*z1(ig) |
---|
| 613 | ENDDO |
---|
| 614 | |
---|
| 615 | ! Convergence loop |
---|
| 616 | |
---|
| 617 | DO j=1,llnt(ig) |
---|
| 618 | |
---|
| 619 | z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) |
---|
| 620 | zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay) |
---|
| 621 | & -betam(ig,nlay)*dmice(ig,nlay) |
---|
| 622 | zc(ig,nlay)=zc(ig,nlay)*z1(ig) |
---|
| 623 | ! zd(ig,nlay)=zb(ig,nlay)*z1(ig) |
---|
| 624 | |
---|
| 625 | ! calcul de zc et zd pour les couches du haut vers le bas |
---|
| 626 | |
---|
| 627 | DO ilay=nlay-1,1,-1 |
---|
| 628 | z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ |
---|
| 629 | $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) |
---|
| 630 | zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ |
---|
| 631 | $ zb(ig,ilay+1)*zc(ig,ilay+1)- |
---|
| 632 | $ betam(ig,ilay)*dmice(ig,ilay))*z1(ig) |
---|
| 633 | ! zd(ig,ilay)=zb(ig,ilay)*z1(ig) |
---|
| 634 | ENDDO |
---|
| 635 | |
---|
[38] | 636 | c ** calcul de la temperature_d'interface et de sa tendance. |
---|
| 637 | c on ecrit que la somme des flux est nulle a l'interface |
---|
| 638 | c a t + \delta t, |
---|
| 639 | c c'est a dire le flux radiatif a {t + \delta t} |
---|
| 640 | c + le flux turbulent a {t + \delta t} |
---|
| 641 | c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1 |
---|
| 642 | c (notation K dt = /cpp*b/) |
---|
| 643 | c + le flux dans le sol a t |
---|
| 644 | c + l'evolution du flux dans le sol lorsque la temperature d'interface |
---|
| 645 | c passe de sa valeur a t a sa valeur a {t + \delta t}. |
---|
| 646 | c ---------------------------------------------------- |
---|
| 647 | |
---|
| 648 | z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) |
---|
| 649 | s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep |
---|
| 650 | z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig) |
---|
| 651 | ztsrf2(ig)=z1(ig)/z2(ig) |
---|
[473] | 652 | ! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop |
---|
| 653 | zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) |
---|
[38] | 654 | |
---|
| 655 | c ** et a partir de la temperature au sol on remonte |
---|
| 656 | c ----------------------------------------------- |
---|
| 657 | |
---|
[473] | 658 | DO ilay=2,nlay |
---|
| 659 | zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1) |
---|
| 660 | ENDDO |
---|
| 661 | DO ilay=1,nlay |
---|
| 662 | zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay) |
---|
| 663 | ENDDO |
---|
| 664 | |
---|
| 665 | c Condensation/sublimation in the atmosphere |
---|
| 666 | c ------------------------------------------ |
---|
| 667 | c (computation of zcondicea and dmice) |
---|
| 668 | |
---|
[1996] | 669 | IF (.NOT. co2clouds) then |
---|
| 670 | DO l=nlay , 1, -1 |
---|
[473] | 671 | IF(zt(ig,l).LT.ztcond(ig,l)) THEN |
---|
| 672 | pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep |
---|
| 673 | zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1)) |
---|
| 674 | & *ccond*pdtc(ig,l) |
---|
| 675 | dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep |
---|
| 676 | END IF |
---|
[1996] | 677 | ENDDO |
---|
| 678 | ELSE |
---|
| 679 | DO l=nlay , 1, -1 |
---|
| 680 | zcondicea(ig,l)= 0.!pcondicea_co2microp(ig,l)* |
---|
| 681 | c & (pplev(ig,l) - pplev(ig,l+1))/g |
---|
| 682 | dmice(ig,l)= 0.!dmice(ig,l) + zcondicea(ig,l)*ptimestep |
---|
| 683 | pdtc(ig,l)=0. |
---|
| 684 | ENDDO |
---|
| 685 | ENDIF |
---|
| 686 | |
---|
| 687 | ENDDO!of Do j=1,XXX |
---|
[38] | 688 | |
---|
[1996] | 689 | ENDDO !of Do ig=1,ngrid |
---|
[38] | 690 | |
---|
[473] | 691 | pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep |
---|
[669] | 692 | |
---|
[660] | 693 | DO ig=1,ngrid ! computing sensible heat flux (atm => surface) |
---|
| 694 | sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig)) |
---|
| 695 | ENDDO |
---|
[473] | 696 | |
---|
[38] | 697 | c----------------------------------------------------------------------- |
---|
| 698 | c TRACERS |
---|
| 699 | c ------- |
---|
| 700 | |
---|
| 701 | if(tracer) then |
---|
| 702 | |
---|
| 703 | c Using the wind modified by friction for lifting and sublimation |
---|
| 704 | c ---------------------------------------------------------------- |
---|
| 705 | |
---|
[529] | 706 | ! This is computed above and takes into account surface-atmosphere flux |
---|
| 707 | ! enhancement by subgrid gustiness and atmospheric-stability related |
---|
| 708 | ! variations of transfer coefficients. |
---|
| 709 | |
---|
| 710 | ! DO ig=1,ngrid |
---|
| 711 | ! zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) |
---|
| 712 | ! zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig)) |
---|
| 713 | ! zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig)) |
---|
| 714 | ! ENDDO |
---|
| 715 | |
---|
[38] | 716 | c Calcul du flux vertical au bas de la premiere couche (dust) : |
---|
| 717 | c ----------------------------------------------------------- |
---|
[1047] | 718 | do ig=1,ngrid |
---|
[38] | 719 | rho(ig) = zb0(ig,1) /ptimestep |
---|
| 720 | c zb(ig,1) = 0. |
---|
| 721 | end do |
---|
| 722 | c Dust lifting: |
---|
| 723 | if (lifting) then |
---|
[310] | 724 | #ifndef MESOSCALE |
---|
[38] | 725 | if (doubleq.AND.submicron) then |
---|
| 726 | do ig=1,ngrid |
---|
| 727 | c if(co2ice(ig).lt.1) then |
---|
| 728 | pdqsdif(ig,igcm_dust_mass) = |
---|
| 729 | & -alpha_lift(igcm_dust_mass) |
---|
| 730 | pdqsdif(ig,igcm_dust_number) = |
---|
| 731 | & -alpha_lift(igcm_dust_number) |
---|
| 732 | pdqsdif(ig,igcm_dust_submicron) = |
---|
| 733 | & -alpha_lift(igcm_dust_submicron) |
---|
| 734 | c end if |
---|
| 735 | end do |
---|
| 736 | else if (doubleq) then |
---|
[1974] | 737 | if (dustinjection.eq.0) then !injection scheme 0 (old) |
---|
| 738 | !or 2 (injection in CL) |
---|
| 739 | do ig=1,ngrid |
---|
[1455] | 740 | if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2 |
---|
[38] | 741 | pdqsdif(ig,igcm_dust_mass) = |
---|
| 742 | & -alpha_lift(igcm_dust_mass) |
---|
| 743 | pdqsdif(ig,igcm_dust_number) = |
---|
[520] | 744 | & -alpha_lift(igcm_dust_number) |
---|
| 745 | end if |
---|
[1974] | 746 | end do |
---|
| 747 | elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface |
---|
| 748 | do ig=1,ngrid |
---|
| 749 | if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2 |
---|
[2160] | 750 | IF((ti_injection_sol.LE.local_time(ig)).and. |
---|
| 751 | & (local_time(ig).LE.tf_injection_sol)) THEN |
---|
[1974] | 752 | if (rdstorm) then !Rocket dust storm scheme |
---|
| 753 | pdqsdif(ig,igcm_stormdust_mass) = |
---|
| 754 | & -alpha_lift(igcm_stormdust_mass) |
---|
| 755 | & *dustliftday(ig) |
---|
| 756 | pdqsdif(ig,igcm_stormdust_number) = |
---|
| 757 | & -alpha_lift(igcm_stormdust_number) |
---|
| 758 | & *dustliftday(ig) |
---|
| 759 | pdqsdif(ig,igcm_dust_mass)= 0. |
---|
| 760 | pdqsdif(ig,igcm_dust_number)= 0. |
---|
| 761 | else |
---|
| 762 | pdqsdif(ig,igcm_dust_mass)= |
---|
| 763 | & -dustliftday(ig)* |
---|
| 764 | & alpha_lift(igcm_dust_mass) |
---|
| 765 | pdqsdif(ig,igcm_dust_number)= |
---|
| 766 | & -dustliftday(ig)* |
---|
| 767 | & alpha_lift(igcm_dust_number) |
---|
| 768 | endif |
---|
| 769 | if (submicron) then |
---|
| 770 | pdqsdif(ig,igcm_dust_submicron) = 0. |
---|
| 771 | endif ! if (submicron) |
---|
| 772 | ELSE ! outside dust injection time frame |
---|
[2080] | 773 | pdqsdif(ig,igcm_dust_mass)= 0. |
---|
| 774 | pdqsdif(ig,igcm_dust_number)= 0. |
---|
| 775 | if (rdstorm) then |
---|
[1974] | 776 | pdqsdif(ig,igcm_stormdust_mass)= 0. |
---|
| 777 | pdqsdif(ig,igcm_stormdust_number)= 0. |
---|
[2080] | 778 | end if |
---|
[1974] | 779 | ENDIF |
---|
| 780 | |
---|
| 781 | end if ! of if(co2ice(ig).lt.1) |
---|
| 782 | end do |
---|
| 783 | endif ! end if dustinjection |
---|
[38] | 784 | else if (submicron) then |
---|
| 785 | do ig=1,ngrid |
---|
| 786 | pdqsdif(ig,igcm_dust_submicron) = |
---|
| 787 | & -alpha_lift(igcm_dust_submicron) |
---|
| 788 | end do |
---|
| 789 | else |
---|
[1236] | 790 | #endif |
---|
[38] | 791 | call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice, |
---|
| 792 | & pdqsdif) |
---|
[1236] | 793 | #ifndef MESOSCALE |
---|
[38] | 794 | endif !doubleq.AND.submicron |
---|
[310] | 795 | #endif |
---|
[38] | 796 | else |
---|
| 797 | pdqsdif(1:ngrid,1:nq) = 0. |
---|
| 798 | end if |
---|
| 799 | |
---|
| 800 | c OU calcul de la valeur de q a la surface (water) : |
---|
| 801 | c ---------------------------------------- |
---|
| 802 | |
---|
| 803 | c Inversion pour l'implicite sur q |
---|
[2515] | 804 | c Cas des traceurs qui ne sont pas h2o_vap |
---|
| 805 | c h2o_vap est traite plus loin avec un sous pas de temps |
---|
| 806 | c hdo_vap est traite ensuite car dependant de h2o_vap |
---|
[38] | 807 | c -------------------------------- |
---|
[2515] | 808 | |
---|
| 809 | do iq=1,nq !for all tracers except water vapor |
---|
| 810 | if ((.not. water).or.(.not. iq.eq.igcm_h2o_vap).or. |
---|
| 811 | & (.not. iq.eq.igcm_hdo_vap)) then |
---|
| 812 | |
---|
| 813 | |
---|
[2274] | 814 | zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) |
---|
[2515] | 815 | zb(1:ngrid,1)=0 |
---|
[38] | 816 | |
---|
[2515] | 817 | DO ig=1,ngrid |
---|
| 818 | z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) |
---|
| 819 | zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) |
---|
| 820 | zd(ig,nlay)=zb(ig,nlay)*z1(ig) |
---|
| 821 | ENDDO |
---|
| 822 | |
---|
| 823 | DO ilay=nlay-1,2,-1 |
---|
| 824 | DO ig=1,ngrid |
---|
| 825 | z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ |
---|
| 826 | $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) |
---|
| 827 | zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ |
---|
| 828 | $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) |
---|
| 829 | zd(ig,ilay)=zb(ig,ilay)*z1(ig) |
---|
| 830 | ENDDO |
---|
| 831 | ENDDO |
---|
| 832 | |
---|
| 833 | if ((iq.eq.igcm_h2o_ice) |
---|
| 834 | $ .or. (hdo.and.(iq.eq.igcm_hdo_ice) )) then |
---|
| 835 | |
---|
| 836 | DO ig=1,ngrid |
---|
| 837 | z1(ig)=1./(za(ig,1)+zb(ig,1)+ |
---|
| 838 | $ zb(ig,2)*(1.-zd(ig,2))) |
---|
| 839 | zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ |
---|
| 840 | $ zb(ig,2)*zc(ig,2)) *z1(ig) !special case h2o_ice |
---|
| 841 | ENDDO |
---|
| 842 | else ! every other tracer |
---|
| 843 | DO ig=1,ngrid |
---|
| 844 | z1(ig)=1./(za(ig,1)+zb(ig,1)+ |
---|
| 845 | $ zb(ig,2)*(1.-zd(ig,2))) |
---|
| 846 | zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ |
---|
| 847 | $ zb(ig,2)*zc(ig,2) + |
---|
| 848 | $ (-pdqsdif(ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface |
---|
| 849 | ENDDO |
---|
| 850 | endif !((iq.eq.igcm_h2o_ice) |
---|
| 851 | c Starting upward calculations for simple mixing of tracer (dust) |
---|
| 852 | DO ig=1,ngrid |
---|
| 853 | zq(ig,1,iq)=zc(ig,1) |
---|
| 854 | DO ilay=2,nlay |
---|
| 855 | zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) |
---|
| 856 | ENDDO |
---|
| 857 | ENDDO |
---|
| 858 | endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then |
---|
| 859 | enddo ! of do iq=1,nq |
---|
| 860 | |
---|
| 861 | c --------- h2o_vap -------------------------------- |
---|
| 862 | |
---|
| 863 | |
---|
| 864 | c Traitement de la vapeur d'eau h2o_vap |
---|
| 865 | c Utilisation d'un sous pas de temps afin |
---|
| 866 | c de decrire le flux de chaleur latente |
---|
| 867 | |
---|
| 868 | |
---|
| 869 | do iq=1,nq |
---|
[2312] | 870 | if ((water).and.(iq.eq.igcm_h2o_vap)) then |
---|
[2515] | 871 | |
---|
| 872 | |
---|
| 873 | DO ig=1,ngrid |
---|
| 874 | zqsurf(ig)=pqsurf(ig,igcm_h2o_ice) |
---|
| 875 | ENDDO ! ig=1,ngrid |
---|
| 876 | |
---|
| 877 | c make_tsub : sous pas de temps adaptatif |
---|
| 878 | c la subroutine est a la fin du fichier |
---|
| 879 | |
---|
| 880 | call make_tsub(ngrid,pdtsrf,zqsurf, |
---|
| 881 | & ptimestep,dtmax,watercaptag, |
---|
| 882 | & nsubtimestep) |
---|
| 883 | |
---|
| 884 | c Calculation for turbulent exchange with the surface (for ice) |
---|
| 885 | c initialization of ztsrf, which is surface temperature in |
---|
| 886 | c the subtimestep. |
---|
| 887 | DO ig=1,ngrid |
---|
| 888 | subtimestep = ptimestep/nsubtimestep(ig) |
---|
| 889 | ztsrf(ig)=ptsrf(ig) ! +pdtsrf(ig)*subtimestep |
---|
| 890 | |
---|
| 891 | c Debut du sous pas de temps |
---|
| 892 | |
---|
| 893 | DO tsub=1,nsubtimestep(ig) |
---|
| 894 | |
---|
| 895 | c C'est parti ! |
---|
| 896 | |
---|
| 897 | zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) |
---|
| 898 | & /float(nsubtimestep(ig)) |
---|
[2274] | 899 | zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1) |
---|
[2515] | 900 | & /float(nsubtimestep(ig)) |
---|
[2274] | 901 | zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1) |
---|
[2515] | 902 | |
---|
| 903 | z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) |
---|
| 904 | zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) |
---|
| 905 | zd(ig,nlay)=zb(ig,nlay)*z1(ig) |
---|
| 906 | DO ilay=nlay-1,2,-1 |
---|
| 907 | z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ |
---|
| 908 | $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) |
---|
| 909 | zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ |
---|
| 910 | $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) |
---|
| 911 | zd(ig,ilay)=zb(ig,ilay)*z1(ig) |
---|
| 912 | ENDDO |
---|
| 913 | z1(ig)=1./(za(ig,1)+zb(ig,1)+ |
---|
| 914 | $ zb(ig,2)*(1.-zd(ig,2))) |
---|
| 915 | zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ |
---|
| 916 | $ zb(ig,2)*zc(ig,2)) * z1(ig) |
---|
[38] | 917 | |
---|
[2531] | 918 | call watersat(1,ztsrf(ig),pplev(ig,1),qsat(ig)) |
---|
[2515] | 919 | old_h2o_vap(ig)=zq(ig,1,igcm_h2o_vap) |
---|
| 920 | zd(ig,1)=zb(ig,1)*z1(ig) |
---|
| 921 | zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) |
---|
| 922 | |
---|
| 923 | zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig) |
---|
| 924 | & *(zq1temp(ig)-qsat(ig)) |
---|
[2587] | 925 | c write(*,*)'subliming more than available frost: qsurf!' |
---|
| 926 | if(.not.watercaptag(ig)) then |
---|
| 927 | if ((-zdqsdif(ig)*subtimestep) |
---|
[2515] | 928 | & .gt.(zqsurf(ig))) then |
---|
[2587] | 929 | c pdqsdif > 0 : ice condensing |
---|
| 930 | c pdqsdif < 0 : ice subliming |
---|
| 931 | c write(*,*) "subliming more than available frost: qsurf!" |
---|
[2515] | 932 | zdqsdif(ig)= |
---|
| 933 | & -zqsurf(ig)/subtimestep |
---|
| 934 | c write(*,*)'flux vers le sol=',pdqsdif(ig,nq) |
---|
| 935 | z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2))) |
---|
| 936 | zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+ |
---|
[2587] | 937 | $ zb(ig,2)*zc(ig,2) + |
---|
| 938 | $ (-zdqsdif(ig)) *subtimestep) *z1(ig) |
---|
[2515] | 939 | zq1temp(ig)=zc(ig,1) |
---|
| 940 | endif !if .not.watercaptag(ig) |
---|
| 941 | endif ! if sublim more than surface |
---|
| 942 | |
---|
| 943 | c Starting upward calculations for water : |
---|
| 944 | c Actualisation de h2o_vap dans le premier niveau |
---|
| 945 | zq(ig,1,igcm_h2o_vap)=zq1temp(ig) |
---|
[2587] | 946 | |
---|
| 947 | c Take into account the H2O latent heat impact on the surface temperature |
---|
[2515] | 948 | if (latentheat_surfwater) then |
---|
| 949 | lh=(2834.3-0.28*(ztsrf(ig)-To)- |
---|
| 950 | & 0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3 |
---|
[2587] | 951 | zdtsrf(ig)= zdqsdif(ig)*lh /pcapcal(ig) |
---|
[2515] | 952 | endif ! (latentheat_surfwater) then |
---|
| 953 | |
---|
[2587] | 954 | DO ilay=2,nlay |
---|
| 955 | zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) |
---|
| 956 | ENDDO |
---|
[2515] | 957 | |
---|
| 958 | c Subtimestep water budget : |
---|
| 959 | |
---|
| 960 | ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig) + zdtsrf(ig)) |
---|
| 961 | & *subtimestep |
---|
| 962 | zqsurf(ig)= zqsurf(ig)+( |
---|
| 963 | & zdqsdif(ig))*subtimestep |
---|
| 964 | |
---|
| 965 | |
---|
| 966 | c We ensure that surface temperature can't rise above the solid domain if there |
---|
| 967 | c is still ice on the surface (oldschool) |
---|
| 968 | if(zqsurf(ig) |
---|
| 969 | & +zdqsdif(ig)*subtimestep |
---|
| 970 | & .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To |
---|
| 971 | & zdtsrf(ig) = min(zdtsrf(ig),(To-ztsrf(ig))/subtimestep) ! ice melt case |
---|
| 972 | |
---|
| 973 | |
---|
[2587] | 974 | |
---|
[2515] | 975 | c Fin du sous pas de temps |
---|
| 976 | ENDDO ! tsub=1,nsubtimestep |
---|
| 977 | |
---|
[2587] | 978 | c Integration of subtimestep temp and water budget : |
---|
| 979 | c (btw could also compute the post timestep temp and ice |
---|
| 980 | c by simply adding the subtimestep trend instead of this) |
---|
[2515] | 981 | pdtsrf(ig)= (ztsrf(ig) - |
---|
| 982 | & ptsrf(ig))/ptimestep |
---|
| 983 | pdqsdif(ig,igcm_h2o_ice)= |
---|
| 984 | & (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice))/ptimestep |
---|
| 985 | |
---|
[2587] | 986 | c if subliming more than qsurf(ice) and on watercaptag, water |
---|
| 987 | c sublimates from watercap reservoir (dwatercap_dif is <0) |
---|
| 988 | if(watercaptag(ig)) then |
---|
| 989 | if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep) |
---|
| 990 | & .gt.(pqsurf(ig,igcm_h2o_ice))) then |
---|
| 991 | dwatercap_dif(ig)=pdqsdif(ig,igcm_h2o_ice)+ |
---|
| 992 | & pqsurf(ig,igcm_h2o_ice)/ptimestep |
---|
| 993 | pdqsdif(ig,igcm_h2o_ice)= |
---|
| 994 | & - pqsurf(ig,igcm_h2o_ice)/ptimestep |
---|
| 995 | endif! ((-pdqsdif(ig)*ptimestep) |
---|
| 996 | endif !(watercaptag(ig)) then |
---|
| 997 | |
---|
[2515] | 998 | ENDDO ! of DO ig=1,ngrid |
---|
| 999 | END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap)) |
---|
| 1000 | |
---|
| 1001 | c --------- end of h2o_vap ---------------------------- |
---|
| 1002 | |
---|
| 1003 | c --------- hdo_vap ----------------------------------- |
---|
| 1004 | |
---|
| 1005 | c hdo_ice has already been with along h2o_ice |
---|
| 1006 | c amongst "normal" tracers (ie not h2o_vap) |
---|
| 1007 | |
---|
| 1008 | if (hdo.and.(iq.eq.igcm_hdo_vap)) then |
---|
| 1009 | zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) |
---|
| 1010 | zb(1:ngrid,1)=0 |
---|
| 1011 | |
---|
[38] | 1012 | DO ig=1,ngrid |
---|
| 1013 | z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) |
---|
| 1014 | zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) |
---|
| 1015 | zd(ig,nlay)=zb(ig,nlay)*z1(ig) |
---|
| 1016 | ENDDO |
---|
[2515] | 1017 | |
---|
[38] | 1018 | DO ilay=nlay-1,2,-1 |
---|
| 1019 | DO ig=1,ngrid |
---|
| 1020 | z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ |
---|
| 1021 | $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) |
---|
| 1022 | zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ |
---|
| 1023 | $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) |
---|
| 1024 | zd(ig,ilay)=zb(ig,ilay)*z1(ig) |
---|
| 1025 | ENDDO |
---|
| 1026 | ENDDO |
---|
| 1027 | |
---|
[2312] | 1028 | CALL hdo_surfex(ngrid,nlay,nq,ptimestep, |
---|
[2593] | 1029 | & zt,pplay,zq,pqsurf, |
---|
| 1030 | & old_h2o_vap,qsat,pdqsdif,dwatercap_dif, |
---|
[2312] | 1031 | & hdoflux) |
---|
| 1032 | DO ig=1,ngrid |
---|
| 1033 | z1(ig)=1./(za(ig,1)+zb(ig,1)+ |
---|
| 1034 | $ zb(ig,2)*(1.-zd(ig,2))) |
---|
| 1035 | zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ |
---|
| 1036 | $ zb(ig,2)*zc(ig,2) + |
---|
| 1037 | $ (-hdoflux(ig)) *ptimestep) *z1(ig) !tracer flux from surface |
---|
| 1038 | ENDDO |
---|
| 1039 | |
---|
[38] | 1040 | DO ig=1,ngrid |
---|
[2515] | 1041 | zq(ig,1,iq)=zc(ig,1) |
---|
| 1042 | DO ilay=2,nlay |
---|
| 1043 | zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) |
---|
| 1044 | ENDDO |
---|
[38] | 1045 | ENDDO |
---|
[2515] | 1046 | endif ! (hdo.and.(iq.eq.igcm_hdo_vap)) |
---|
[2312] | 1047 | |
---|
[2515] | 1048 | c --------- end of hdo ---------------------------- |
---|
[38] | 1049 | |
---|
[2515] | 1050 | enddo ! of do iq=1,nq |
---|
| 1051 | end if ! of if(tracer) |
---|
[38] | 1052 | |
---|
[2515] | 1053 | c --------- end of tracers ---------------------------- |
---|
[2312] | 1054 | |
---|
[2260] | 1055 | |
---|
[2312] | 1056 | C Diagnostic output for HDO |
---|
| 1057 | if (hdo) then |
---|
| 1058 | CALL WRITEDIAGFI(ngrid,'hdoflux', |
---|
| 1059 | & 'hdoflux', |
---|
| 1060 | & ' ',2,hdoflux) |
---|
| 1061 | CALL WRITEDIAGFI(ngrid,'h2oflux', |
---|
| 1062 | & 'h2oflux', |
---|
| 1063 | & ' ',2,h2oflux) |
---|
| 1064 | endif |
---|
| 1065 | |
---|
[38] | 1066 | c----------------------------------------------------------------------- |
---|
| 1067 | c 8. calcul final des tendances de la diffusion verticale |
---|
| 1068 | c ---------------------------------------------------- |
---|
| 1069 | |
---|
| 1070 | DO ilev = 1, nlay |
---|
| 1071 | DO ig=1,ngrid |
---|
| 1072 | pdudif(ig,ilev)=( zu(ig,ilev)- |
---|
| 1073 | $ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep |
---|
| 1074 | pdvdif(ig,ilev)=( zv(ig,ilev)- |
---|
| 1075 | $ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep |
---|
[473] | 1076 | hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep |
---|
| 1077 | $ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev) |
---|
| 1078 | pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep |
---|
[38] | 1079 | ENDDO |
---|
| 1080 | ENDDO |
---|
| 1081 | |
---|
| 1082 | if (tracer) then |
---|
| 1083 | DO iq = 1, nq |
---|
| 1084 | DO ilev = 1, nlay |
---|
| 1085 | DO ig=1,ngrid |
---|
| 1086 | pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)- |
---|
| 1087 | $ (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep |
---|
| 1088 | ENDDO |
---|
| 1089 | ENDDO |
---|
| 1090 | ENDDO |
---|
| 1091 | end if |
---|
| 1092 | |
---|
| 1093 | c ** diagnostique final |
---|
| 1094 | c ------------------ |
---|
| 1095 | |
---|
| 1096 | IF(lecrit) THEN |
---|
| 1097 | PRINT*,'In vdif' |
---|
| 1098 | PRINT*,'Ts (t) and Ts (t+st)' |
---|
| 1099 | WRITE(*,'(a10,3a15)') |
---|
| 1100 | s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)' |
---|
| 1101 | PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1) |
---|
| 1102 | DO ilev=1,nlay |
---|
| 1103 | WRITE(*,'(4f15.7)') |
---|
[473] | 1104 | s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev), |
---|
[38] | 1105 | s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev) |
---|
| 1106 | |
---|
| 1107 | ENDDO |
---|
| 1108 | ENDIF |
---|
| 1109 | |
---|
[1036] | 1110 | END SUBROUTINE vdifc |
---|
[1969] | 1111 | |
---|
[2312] | 1112 | c==================================== |
---|
| 1113 | |
---|
[2515] | 1114 | SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep, |
---|
| 1115 | $ dtmax,watercaptag,ntsub) |
---|
[2312] | 1116 | |
---|
[2515] | 1117 | c Pas de temps adaptatif en estimant le taux de sublimation |
---|
| 1118 | c et en adaptant avec un critere "dtmax" du chauffage a accomoder |
---|
| 1119 | c dtmax est regle empiriquement (pour l'instant) a 0.5 K |
---|
| 1120 | |
---|
| 1121 | integer,intent(in) :: naersize |
---|
| 1122 | real,intent(in) :: dtsurf(naersize) |
---|
| 1123 | real,intent(in) :: qsurf(naersize) |
---|
| 1124 | logical,intent(in) :: watercaptag(naersize) |
---|
| 1125 | real,intent(in) :: ptimestep |
---|
| 1126 | real,intent(in) :: dtmax |
---|
| 1127 | real :: ztsub(naersize) |
---|
| 1128 | integer :: i |
---|
| 1129 | integer,intent(out) :: ntsub(naersize) |
---|
| 1130 | |
---|
| 1131 | do i=1,naersize |
---|
| 1132 | if ((qsurf(i).eq.0).and. |
---|
| 1133 | & (.not.watercaptag(i))) then |
---|
| 1134 | ntsub(i) = 1 |
---|
| 1135 | else |
---|
| 1136 | ztsub(i) = ptimestep * dtsurf(i) / dtmax |
---|
| 1137 | ntsub(i) = ceiling(abs(ztsub(i))) |
---|
| 1138 | endif ! (qsurf(i).eq.0) then |
---|
| 1139 | c |
---|
| 1140 | c write(78,*), dtsurf*ptimestep, dtsurf, ntsub |
---|
| 1141 | enddo! 1=1,ngrid |
---|
| 1142 | |
---|
| 1143 | |
---|
| 1144 | |
---|
| 1145 | END SUBROUTINE make_tsub |
---|
[1969] | 1146 | END MODULE vdifc_mod |
---|