| 1 | subroutine condense_n2(klon,klev,nq,ptimestep, & |
|---|
| 2 | pcapcal,pplay,pplev,ptsrf,pt, & |
|---|
| 3 | pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, & |
|---|
| 4 | picen2,psolaralb,pemisurf, & |
|---|
| 5 | pdtc,pdtsrfc,pdpsrf,pduc,pdvc, & |
|---|
| 6 | pdqc,pdicen2) |
|---|
| 7 | |
|---|
| 8 | use radinc_h, only : naerkind |
|---|
| 9 | use comgeomfi_h |
|---|
| 10 | use comcstfi_mod, only: g, r, cpp, pi |
|---|
| 11 | USE surfdat_h, only: phisfi,kp,p00 |
|---|
| 12 | USE tracer_h, only: noms, igcm_n2, lw_n2 |
|---|
| 13 | USE callkeys_mod, only: fast,ch4lag,latlag,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag |
|---|
| 14 | USE comvert_mod, ONLY: ap,bp |
|---|
| 15 | use geometry_mod, only: latitude |
|---|
| 16 | |
|---|
| 17 | |
|---|
| 18 | implicit none |
|---|
| 19 | |
|---|
| 20 | !================================================================== |
|---|
| 21 | ! Purpose |
|---|
| 22 | ! ------- |
|---|
| 23 | ! Condense and/or sublime N2 ice on the ground and in the |
|---|
| 24 | ! atmosphere, and sediment the ice. |
|---|
| 25 | ! |
|---|
| 26 | ! Inputs |
|---|
| 27 | ! ------ |
|---|
| 28 | ! klon Number of vertical columns |
|---|
| 29 | ! klev Number of layers |
|---|
| 30 | ! pplay(klon,klev) Pressure layers |
|---|
| 31 | ! pplev(klon,klev+1) Pressure levels |
|---|
| 32 | ! pt(klon,klev) Temperature (in K) |
|---|
| 33 | ! ptsrf(klon) Surface temperature |
|---|
| 34 | ! |
|---|
| 35 | ! pdt(klon,klev) Time derivative before condensation/sublimation of pt |
|---|
| 36 | ! pdtsrf(klon) Time derivative before condensation/sublimation of ptsrf |
|---|
| 37 | ! picen2(klon) n2 ice at the surface (kg/m2) |
|---|
| 38 | ! |
|---|
| 39 | ! Outputs |
|---|
| 40 | ! ------- |
|---|
| 41 | ! pdpsrf(klon) \ Contribution of condensation/sublimation |
|---|
| 42 | ! pdtc(klon,klev) / to the time derivatives of Ps, pt, and ptsrf |
|---|
| 43 | ! pdtsrfc(klon) / |
|---|
| 44 | ! pdicen2(klon) Tendancy n2 ice at the surface (kg/m2) |
|---|
| 45 | ! |
|---|
| 46 | ! Both |
|---|
| 47 | ! ---- |
|---|
| 48 | ! psolaralb(klon) Albedo at the surface |
|---|
| 49 | ! pemisurf(klon) Emissivity of the surface |
|---|
| 50 | ! |
|---|
| 51 | ! Authors |
|---|
| 52 | ! ------- |
|---|
| 53 | ! Francois Forget (1996,2013) |
|---|
| 54 | ! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009) |
|---|
| 55 | ! |
|---|
| 56 | ! |
|---|
| 57 | !================================================================== |
|---|
| 58 | |
|---|
| 59 | !----------------------------------------------------------------------- |
|---|
| 60 | ! Arguments |
|---|
| 61 | |
|---|
| 62 | INTEGER klon, klev, nq |
|---|
| 63 | |
|---|
| 64 | REAL ptimestep |
|---|
| 65 | REAL pplay(klon,klev),pplev(klon,klev+1) |
|---|
| 66 | REAL pcapcal(klon) |
|---|
| 67 | REAL pt(klon,klev) |
|---|
| 68 | REAL ptsrf(klon),flu1(klon),flu2(klon),flu3(klon) |
|---|
| 69 | REAL pphi(klon,klev) |
|---|
| 70 | REAL pdt(klon,klev),pdtsrf(klon),pdtc(klon,klev) |
|---|
| 71 | REAL pdtsrfc(klon),pdpsrf(klon) |
|---|
| 72 | REAL picen2(klon),psolaralb(klon),pemisurf(klon) |
|---|
| 73 | |
|---|
| 74 | |
|---|
| 75 | REAL pu(klon,klev) , pv(klon,klev) |
|---|
| 76 | REAL pdu(klon,klev) , pdv(klon,klev) |
|---|
| 77 | REAL pduc(klon,klev) , pdvc(klon,klev) |
|---|
| 78 | REAL pq(klon,klev,nq),pdq(klon,klev,nq) |
|---|
| 79 | REAL pdqc(klon,klev,nq) |
|---|
| 80 | |
|---|
| 81 | !----------------------------------------------------------------------- |
|---|
| 82 | ! Local variables |
|---|
| 83 | |
|---|
| 84 | INTEGER l,ig,ilay,it,iq,ich4_gas |
|---|
| 85 | |
|---|
| 86 | REAL*8 zt(klon,klev) |
|---|
| 87 | REAL tcond_n2 |
|---|
| 88 | REAL pcond_n2 |
|---|
| 89 | REAL glob_average2d ! function |
|---|
| 90 | REAL zqn2(klon,klev) ! N2 MMR used to compute Tcond/zqn2 |
|---|
| 91 | REAL ztcond (klon,klev) |
|---|
| 92 | REAL ztcondsol(klon),zfallheat |
|---|
| 93 | REAL pdicen2(klon) |
|---|
| 94 | REAL zcondicea(klon,klev), zcondices(klon) |
|---|
| 95 | REAL zfallice(klon,klev+1) |
|---|
| 96 | REAL zmflux(klev+1) |
|---|
| 97 | REAL zu(klev),zv(klev) |
|---|
| 98 | REAL zq(klev,nq),zq1(klev) |
|---|
| 99 | REAL ztsrf(klon) |
|---|
| 100 | REAL ztc(klev), ztm(klev+1) |
|---|
| 101 | REAL zum(klev+1) , zvm(klev+1) |
|---|
| 102 | REAL zqm(klev+1,nq),zqm1(klev+1) |
|---|
| 103 | LOGICAL condsub(klon) |
|---|
| 104 | REAL subptimestep |
|---|
| 105 | Integer Ntime |
|---|
| 106 | real masse (klev),w(klev+1) |
|---|
| 107 | real wq(klon,klev+1) |
|---|
| 108 | real vstokes,reff |
|---|
| 109 | real dWtotsn2 |
|---|
| 110 | real condnconsn2(klon) |
|---|
| 111 | real nconsMAXn2 |
|---|
| 112 | ! Special diagnostic variables |
|---|
| 113 | real tconda1(klon,klev) |
|---|
| 114 | real tconda2(klon,klev) |
|---|
| 115 | real zdtsig (klon,klev) |
|---|
| 116 | real zdtlatent (klon,klev) |
|---|
| 117 | real zdt (klon,klev) |
|---|
| 118 | REAL albediceF(klon) |
|---|
| 119 | ! SAVE albediceF |
|---|
| 120 | INTEGER nsubtimestep,itsub !number of subtimestep when calling vl1d |
|---|
| 121 | REAL subtimestep !ptimestep/nsubtimestep |
|---|
| 122 | REAL dtmax |
|---|
| 123 | |
|---|
| 124 | REAL zplevhist(klon) |
|---|
| 125 | REAL zplevnew(klon) |
|---|
| 126 | REAL zplev(klon) |
|---|
| 127 | REAL zpicen2(klon) |
|---|
| 128 | REAL ztsrfhist(klon) |
|---|
| 129 | REAL zdtsrf(klon) |
|---|
| 130 | REAL globzplevnew |
|---|
| 131 | |
|---|
| 132 | REAL vmrn2(klon) |
|---|
| 133 | ! SAVE vmrn2 |
|---|
| 134 | REAL stephan |
|---|
| 135 | DATA stephan/5.67e-08/ ! Stephan Boltzman constant |
|---|
| 136 | SAVE stephan |
|---|
| 137 | !----------------------------------------------------------------------- |
|---|
| 138 | ! Saved local variables |
|---|
| 139 | |
|---|
| 140 | ! REAL latcond |
|---|
| 141 | REAL ccond |
|---|
| 142 | REAL cpice ! for atm condensation |
|---|
| 143 | SAVE cpice |
|---|
| 144 | ! SAVE latcond,ccond |
|---|
| 145 | SAVE ccond |
|---|
| 146 | |
|---|
| 147 | LOGICAL firstcall |
|---|
| 148 | SAVE firstcall |
|---|
| 149 | REAL SSUM |
|---|
| 150 | EXTERNAL SSUM |
|---|
| 151 | |
|---|
| 152 | ! DATA latcond /2.5e5/ |
|---|
| 153 | ! DATA latcond /1.98e5/ |
|---|
| 154 | DATA cpice /1300./ |
|---|
| 155 | DATA firstcall/.true./ |
|---|
| 156 | |
|---|
| 157 | INTEGER,SAVE :: i_n2ice=0 ! n2 ice |
|---|
| 158 | CHARACTER(LEN=20) :: tracername ! to temporarily store text |
|---|
| 159 | logical olkin ! option to prevent N2 ice effect in the south |
|---|
| 160 | DATA olkin/.false./ |
|---|
| 161 | save olkin |
|---|
| 162 | |
|---|
| 163 | CHARACTER(len=10) :: tname |
|---|
| 164 | |
|---|
| 165 | !----------------------------------------------------------------------- |
|---|
| 166 | |
|---|
| 167 | ! Initialisation |
|---|
| 168 | IF (firstcall) THEN |
|---|
| 169 | ccond=cpp/(g*lw_n2) |
|---|
| 170 | print*,'In condense_n2cloud: ccond=',ccond,' latcond=',lw_n2 |
|---|
| 171 | |
|---|
| 172 | ! calculate global mean surface pressure for the fast mode |
|---|
| 173 | IF (fast) THEN |
|---|
| 174 | IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon)) |
|---|
| 175 | DO ig=1,klon |
|---|
| 176 | kp(ig) = exp(-phisfi(ig)/(r*38.)) |
|---|
| 177 | ENDDO |
|---|
| 178 | p00=glob_average2d(kp) ! mean pres at ref level |
|---|
| 179 | ENDIF |
|---|
| 180 | |
|---|
| 181 | vmrn2(:) = 1. |
|---|
| 182 | !IF (ch4lag) then |
|---|
| 183 | ! DO ig=1,klon |
|---|
| 184 | ! if (latitude(ig)*180./pi.ge.latlag) then |
|---|
| 185 | ! vmrn2(ig) = vmrlag |
|---|
| 186 | ! endif |
|---|
| 187 | ! ENDDO |
|---|
| 188 | !ENDIF |
|---|
| 189 | !IF (no_n2frost) then |
|---|
| 190 | ! DO ig=1,klon |
|---|
| 191 | ! if (picen2(ig).eq.0.) then |
|---|
| 192 | ! vmrn2(ig) = 1.e-15 |
|---|
| 193 | ! endif |
|---|
| 194 | ! ENDDO |
|---|
| 195 | !ENDIF |
|---|
| 196 | firstcall=.false. |
|---|
| 197 | ENDIF |
|---|
| 198 | |
|---|
| 199 | !----------------------------------------------------------------------- |
|---|
| 200 | ! Calculation of n2 condensation / sublimation |
|---|
| 201 | |
|---|
| 202 | ! Variables used: |
|---|
| 203 | ! picen2(klon) : amount of n2 ice on the ground (kg/m2) |
|---|
| 204 | ! zcondicea(klon,klev): condensation rate in layer l (kg/m2/s) |
|---|
| 205 | ! zcondices(klon) : condensation rate on the ground (kg/m2/s) |
|---|
| 206 | ! zfallice(klon,klev) : amount of ice falling from layer l (kg/m2/s) |
|---|
| 207 | ! zdtlatent(klon,klev): dT/dt due to phase changes (K/s) |
|---|
| 208 | |
|---|
| 209 | ! Tendencies initially set to 0 |
|---|
| 210 | zcondices(1:klon) = 0. |
|---|
| 211 | pdtsrfc(1:klon) = 0. |
|---|
| 212 | pdpsrf(1:klon) = 0. |
|---|
| 213 | ztsrfhist(1:klon) = 0. |
|---|
| 214 | condsub(1:klon) = .false. |
|---|
| 215 | pdicen2(1:klon) = 0. |
|---|
| 216 | zfallheat=0 |
|---|
| 217 | pdqc(1:klon,1:klev,1:nq)=0 |
|---|
| 218 | pdtc(1:klon,1:klev)=0 |
|---|
| 219 | pduc(1:klon,1:klev)=0 |
|---|
| 220 | pdvc(1:klon,1:klev)=0 |
|---|
| 221 | zfallice(1:klon,1:klev+1)=0 |
|---|
| 222 | zcondicea(1:klon,1:klev)=0 |
|---|
| 223 | zdtlatent(1:klon,1:klev)=0 |
|---|
| 224 | zt(1:klon,1:klev)=0. |
|---|
| 225 | |
|---|
| 226 | !----------------------------------------------------------------------- |
|---|
| 227 | ! Atmospheric condensation |
|---|
| 228 | |
|---|
| 229 | ! Condensation / sublimation in the atmosphere |
|---|
| 230 | ! -------------------------------------------- |
|---|
| 231 | ! (calcul of zcondicea , zfallice and pdtc) |
|---|
| 232 | |
|---|
| 233 | zt(1:klon,1:klev)=pt(1:klon,1:klev)+ pdt(1:klon,1:klev)*ptimestep |
|---|
| 234 | if (igcm_n2.ne.0) then |
|---|
| 235 | zqn2(1:klon,1:klev) = 1. ! & temporaire |
|---|
| 236 | ! zqn2(1:klon,1:klev)= pq(1:klon,1:klev,igcm_n2) + pdq(1:klon,1:klev,igcm_n2)*ptimestep |
|---|
| 237 | else |
|---|
| 238 | zqn2(1:klon,1:klev) = 1. |
|---|
| 239 | end if |
|---|
| 240 | |
|---|
| 241 | if (.not.fast) then |
|---|
| 242 | ! Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2 |
|---|
| 243 | DO l=1,klev |
|---|
| 244 | DO ig=1,klon |
|---|
| 245 | ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l)) |
|---|
| 246 | ENDDO |
|---|
| 247 | ENDDO |
|---|
| 248 | |
|---|
| 249 | DO l=klev,1,-1 |
|---|
| 250 | DO ig=1,klon |
|---|
| 251 | pdtc(ig,l)=0. ! final tendancy temperature set to 0 |
|---|
| 252 | |
|---|
| 253 | IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN |
|---|
| 254 | condsub(ig)=.true. !condensation at level l |
|---|
| 255 | IF (zfallice(ig,l+1).gt.0) then |
|---|
| 256 | zfallheat=zfallice(ig,l+1)*& |
|---|
| 257 | (pphi(ig,l+1)-pphi(ig,l) +& |
|---|
| 258 | cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2 |
|---|
| 259 | ELSE |
|---|
| 260 | zfallheat=0. |
|---|
| 261 | ENDIF |
|---|
| 262 | zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep |
|---|
| 263 | zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))& |
|---|
| 264 | *ccond*zdtlatent(ig,l)- zfallheat |
|---|
| 265 | ! Case when the ice from above sublimes entirely |
|---|
| 266 | ! """"""""""""""""""""""""""""""""""""""""""""""" |
|---|
| 267 | IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) & |
|---|
| 268 | .AND. (zfallice(ig,l+1).gt.0)) THEN |
|---|
| 269 | |
|---|
| 270 | zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/& |
|---|
| 271 | (ccond*(pplev(ig,l)-pplev(ig,l+1))) |
|---|
| 272 | zcondicea(ig,l)= -zfallice(ig,l+1) |
|---|
| 273 | END IF |
|---|
| 274 | |
|---|
| 275 | zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1) |
|---|
| 276 | |
|---|
| 277 | END IF |
|---|
| 278 | |
|---|
| 279 | ENDDO |
|---|
| 280 | ENDDO |
|---|
| 281 | endif ! not fast |
|---|
| 282 | |
|---|
| 283 | !----------------------------------------------------------------------- |
|---|
| 284 | ! Condensation/sublimation on the ground |
|---|
| 285 | ! (calculation of zcondices and pdtsrfc) |
|---|
| 286 | |
|---|
| 287 | ! Adding subtimesteps : in fast version, pressures too low lead to |
|---|
| 288 | ! instabilities. |
|---|
| 289 | IF (fast) THEN |
|---|
| 290 | IF (pplev(1,1).gt.0.3) THEN |
|---|
| 291 | nsubtimestep= 1 |
|---|
| 292 | ELSE |
|---|
| 293 | nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1) |
|---|
| 294 | ENDIF |
|---|
| 295 | ELSE |
|---|
| 296 | nsubtimestep= 1 ! if more then kp and p00 have to be calculated |
|---|
| 297 | ! for nofast mode |
|---|
| 298 | ENDIF |
|---|
| 299 | subtimestep=ptimestep/float(nsubtimestep) |
|---|
| 300 | |
|---|
| 301 | do itsub=1,nsubtimestep |
|---|
| 302 | ! first loop : getting zplev, ztsurf |
|---|
| 303 | IF (itsub.eq.1) then |
|---|
| 304 | DO ig=1,klon |
|---|
| 305 | zplev(ig)=pplev(ig,1) |
|---|
| 306 | ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep |
|---|
| 307 | ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep !! |
|---|
| 308 | zpicen2(ig)=picen2(ig) |
|---|
| 309 | ENDDO |
|---|
| 310 | ELSE |
|---|
| 311 | ! additional loop : |
|---|
| 312 | ! 1) pressure updated |
|---|
| 313 | ! 2) direct redistribution of pressure over the globe |
|---|
| 314 | ! 3) modification pressure for unstable cases |
|---|
| 315 | ! 4) pressure update to conserve mass |
|---|
| 316 | ! 5) temperature updated with radiative tendancies |
|---|
| 317 | DO ig=1,klon |
|---|
| 318 | zplevhist(ig)=zplev(ig) |
|---|
| 319 | zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep ! 1) |
|---|
| 320 | !IF (zplevnew(ig).le.0.0001) then |
|---|
| 321 | ! zplevnew(ig)=0.0001*kp(ig)/p00 |
|---|
| 322 | !ENDIF |
|---|
| 323 | ENDDO |
|---|
| 324 | ! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code) |
|---|
| 325 | globzplevnew=glob_average2d(zplevnew) |
|---|
| 326 | DO ig=1,klon |
|---|
| 327 | zplev(ig)=kp(ig)*globzplevnew/p00 ! 2) |
|---|
| 328 | ENDDO |
|---|
| 329 | DO ig=1,klon ! 3) unstable case condensation and sublimation: zplev=zplevhist |
|---|
| 330 | IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or. & |
|---|
| 331 | ((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then |
|---|
| 332 | zplev(ig)=zplevhist(ig) |
|---|
| 333 | ENDIF |
|---|
| 334 | zplevhist(ig)=zplev(ig) |
|---|
| 335 | ENDDO |
|---|
| 336 | zplev=zplev*globzplevnew/glob_average2d(zplevhist) ! 4) |
|---|
| 337 | DO ig=1,klon ! 5) |
|---|
| 338 | zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))*(ptsrf(ig)**4-ztsrf(ig)**4) |
|---|
| 339 | ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)*subtimestep+zdtsrf(ig)*subtimestep |
|---|
| 340 | zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep |
|---|
| 341 | ENDDO |
|---|
| 342 | ENDIF ! (itsub=1) |
|---|
| 343 | |
|---|
| 344 | DO ig=1,klon |
|---|
| 345 | ! forecast of frost temperature ztcondsol |
|---|
| 346 | ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1)) |
|---|
| 347 | !ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig)) |
|---|
| 348 | |
|---|
| 349 | ! Loop over where we have condensation / sublimation |
|---|
| 350 | IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground cond |
|---|
| 351 | ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublim |
|---|
| 352 | (zpicen2(ig) .GT. 0.))) THEN |
|---|
| 353 | condsub(ig) = .true. ! condensation or sublimation |
|---|
| 354 | |
|---|
| 355 | ! Condensation or partial sublimation of N2 ice |
|---|
| 356 | if (ztsrf(ig) .LT. ztcondsol(ig)) then ! condensation |
|---|
| 357 | ! Include a correction to account for the cooling of air near |
|---|
| 358 | ! the surface before condensing: |
|---|
| 359 | zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) & |
|---|
| 360 | /((lw_n2+cpp*(zt(ig,1)-ztcondsol(ig)))*subtimestep) |
|---|
| 361 | else ! sublimation |
|---|
| 362 | zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) & |
|---|
| 363 | /(lw_n2*subtimestep) |
|---|
| 364 | end if |
|---|
| 365 | |
|---|
| 366 | pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep |
|---|
| 367 | |
|---|
| 368 | ! partial sublimation of N2 ice |
|---|
| 369 | ! If the entire N_2 ice layer sublimes |
|---|
| 370 | ! (including what has just condensed in the atmosphere) |
|---|
| 371 | IF((zpicen2(ig)/subtimestep).LE. & |
|---|
| 372 | -zcondices(ig))THEN |
|---|
| 373 | zcondices(ig) = -zpicen2(ig)/subtimestep |
|---|
| 374 | pdtsrfc(ig)=(lw_n2/pcapcal(ig))* & |
|---|
| 375 | (zcondices(ig)) |
|---|
| 376 | END IF |
|---|
| 377 | |
|---|
| 378 | ! Changing N2 ice amount and pressure |
|---|
| 379 | |
|---|
| 380 | pdicen2(ig) = zcondices(ig) |
|---|
| 381 | pdpsrf(ig) = -pdicen2(ig)*g |
|---|
| 382 | ! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond |
|---|
| 383 | IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then |
|---|
| 384 | pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep |
|---|
| 385 | pdicen2(ig)=-pdpsrf(ig)/g |
|---|
| 386 | ENDIF |
|---|
| 387 | |
|---|
| 388 | ELSE ! no condsub |
|---|
| 389 | pdpsrf(ig)=0. |
|---|
| 390 | pdicen2(ig)=0. |
|---|
| 391 | pdtsrfc(ig)=0. |
|---|
| 392 | ENDIF |
|---|
| 393 | ENDDO ! ig |
|---|
| 394 | enddo ! subtimestep |
|---|
| 395 | |
|---|
| 396 | ! Updating pressure, temperature and ice reservoir |
|---|
| 397 | DO ig=1,klon |
|---|
| 398 | pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep |
|---|
| 399 | ! Two options here : 1 ok, 2 is wrong |
|---|
| 400 | pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep |
|---|
| 401 | !pdicen2(ig)=-pdpsrf(ig)/g |
|---|
| 402 | |
|---|
| 403 | pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep |
|---|
| 404 | |
|---|
| 405 | ! security |
|---|
| 406 | if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then |
|---|
| 407 | write(*,*) 'WARNING in condense_n2:' |
|---|
| 408 | write(*,*) picen2(ig),pdicen2(ig)*ptimestep |
|---|
| 409 | pdicen2(ig)= -picen2(ig)/ptimestep |
|---|
| 410 | pdpsrf(ig)=-pdicen2(ig)*g |
|---|
| 411 | endif |
|---|
| 412 | |
|---|
| 413 | if(.not.picen2(ig).ge.0.) THEN |
|---|
| 414 | ! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then |
|---|
| 415 | print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep |
|---|
| 416 | ! pdicen2(ig)= -picen2(ig)/ptimestep |
|---|
| 417 | ! else |
|---|
| 418 | picen2(ig)=0.0 |
|---|
| 419 | ! endif |
|---|
| 420 | endif |
|---|
| 421 | ENDDO |
|---|
| 422 | |
|---|
| 423 | ! *************************************************************** |
|---|
| 424 | ! Correction to account for redistribution between sigma or hybrid |
|---|
| 425 | ! layers when changing surface pressure (and warming/cooling |
|---|
| 426 | ! of the n2 currently changing phase). |
|---|
| 427 | ! ************************************************************* |
|---|
| 428 | if (.not.fast) then |
|---|
| 429 | DO ig=1,klon |
|---|
| 430 | if (condsub(ig)) then |
|---|
| 431 | |
|---|
| 432 | ! Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) |
|---|
| 433 | ! """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" |
|---|
| 434 | zmflux(1) = -zcondices(ig) |
|---|
| 435 | DO l=1,klev |
|---|
| 436 | zmflux(l+1) = zmflux(l) -zcondicea(ig,l) & |
|---|
| 437 | + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1)) |
|---|
| 438 | ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld |
|---|
| 439 | if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0. |
|---|
| 440 | END DO |
|---|
| 441 | |
|---|
| 442 | ! Mass of each layer |
|---|
| 443 | ! ------------------ |
|---|
| 444 | DO l=1,klev |
|---|
| 445 | masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g |
|---|
| 446 | END DO |
|---|
| 447 | |
|---|
| 448 | |
|---|
| 449 | ! Corresponding fluxes for T,U,V,Q |
|---|
| 450 | ! """""""""""""""""""""""""""""""" |
|---|
| 451 | ! averaging operator for TRANSPORT |
|---|
| 452 | ! """""""""""""""""""""""""""""""" |
|---|
| 453 | |
|---|
| 454 | ! Subtimestep loop to perform the redistribution gently and simultaneously with |
|---|
| 455 | ! the other tendencies |
|---|
| 456 | ! Estimation of subtimestep (using only the first layer, the most critical) |
|---|
| 457 | dtmax=ptimestep |
|---|
| 458 | if (zmflux(1).gt.1.e-20) then |
|---|
| 459 | dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1))) |
|---|
| 460 | endif |
|---|
| 461 | nsubtimestep= max(nint(ptimestep/dtmax),nint(2.)) |
|---|
| 462 | subtimestep=ptimestep/float(nsubtimestep) |
|---|
| 463 | |
|---|
| 464 | ! New flux for each subtimestep |
|---|
| 465 | do l=1,klev+1 |
|---|
| 466 | w(l)=-zmflux(l)*subtimestep |
|---|
| 467 | enddo |
|---|
| 468 | ! initializing variables that will vary during subtimestep: |
|---|
| 469 | do l=1,klev |
|---|
| 470 | ztc(l) =pt(ig,l) |
|---|
| 471 | zu(l) =pu(ig,l) |
|---|
| 472 | zv(l) =pv(ig,l) |
|---|
| 473 | do iq=1,nq |
|---|
| 474 | zq(l,iq) = pq(ig,l,iq) |
|---|
| 475 | enddo |
|---|
| 476 | end do |
|---|
| 477 | |
|---|
| 478 | ! loop over nsubtimestep |
|---|
| 479 | ! """""""""""""""""""""" |
|---|
| 480 | do itsub=1,nsubtimestep |
|---|
| 481 | ! Progressively adding tendancies from other processes. |
|---|
| 482 | do l=1,klev |
|---|
| 483 | ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep |
|---|
| 484 | zu(l) =zu(l) +pdu( ig,l) * subtimestep |
|---|
| 485 | zv(l) =zv(l) +pdv( ig,l) * subtimestep |
|---|
| 486 | do iq=1,nq |
|---|
| 487 | zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep |
|---|
| 488 | enddo |
|---|
| 489 | end do |
|---|
| 490 | |
|---|
| 491 | ! Change of mass in each layer |
|---|
| 492 | do l=1,klev |
|---|
| 493 | masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))& |
|---|
| 494 | /(g*pplev(ig,1)) |
|---|
| 495 | end do |
|---|
| 496 | |
|---|
| 497 | ! Value transfert at the surface interface when condensation sublimation: |
|---|
| 498 | |
|---|
| 499 | if (zmflux(1).lt.0) then |
|---|
| 500 | ! Surface condensation |
|---|
| 501 | zum(1)= zu(1) |
|---|
| 502 | zvm(1)= zv(1) |
|---|
| 503 | ztm(1) = ztc(1) |
|---|
| 504 | else |
|---|
| 505 | ! Surface sublimation: |
|---|
| 506 | ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep |
|---|
| 507 | zum(1) = 0 |
|---|
| 508 | zvm(1) = 0 |
|---|
| 509 | end if |
|---|
| 510 | do iq=1,nq |
|---|
| 511 | zqm(1,iq)=0. ! most tracer do not condense ! |
|---|
| 512 | enddo |
|---|
| 513 | ! Special case if the tracer is n2 gas |
|---|
| 514 | if (igcm_n2.ne.0) zqm(1,igcm_n2)=1. |
|---|
| 515 | |
|---|
| 516 | ztm(2:klev+1)=0. |
|---|
| 517 | zum(2:klev+1)=0. |
|---|
| 518 | zvm(2:klev+1)=0. |
|---|
| 519 | zqm1(1:klev+1)=0. |
|---|
| 520 | |
|---|
| 521 | ! Van Leer scheme: |
|---|
| 522 | call vl1d(klev,ztc,2.,masse,w,ztm) |
|---|
| 523 | call vl1d(klev,zu ,2.,masse,w,zum) |
|---|
| 524 | call vl1d(klev,zv ,2.,masse,w,zvm) |
|---|
| 525 | do iq=1,nq |
|---|
| 526 | do l=1,klev |
|---|
| 527 | zq1(l)=zq(l,iq) |
|---|
| 528 | enddo |
|---|
| 529 | zqm1(1)=zqm(1,iq) |
|---|
| 530 | call vl1d(klev,zq1,2.,masse,w,zqm1) |
|---|
| 531 | do l=2,klev |
|---|
| 532 | zqm(l,iq)=zqm1(l) |
|---|
| 533 | enddo |
|---|
| 534 | enddo |
|---|
| 535 | |
|---|
| 536 | ! Correction to prevent negative value for qn2 |
|---|
| 537 | if (igcm_n2.ne.0) then |
|---|
| 538 | zqm(1,igcm_n2)=1. |
|---|
| 539 | do l=1,klev-1 |
|---|
| 540 | if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then |
|---|
| 541 | zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), & |
|---|
| 542 | (zqm(l,igcm_n2)*w(l) -zq(l,igcm_n2)*masse(l))/w(l+1) ) |
|---|
| 543 | else |
|---|
| 544 | exit |
|---|
| 545 | endif |
|---|
| 546 | end do |
|---|
| 547 | end if |
|---|
| 548 | |
|---|
| 549 | ! Value transfert at the surface interface when condensation sublimation: |
|---|
| 550 | if (zmflux(1).lt.0) then |
|---|
| 551 | ! Surface condensation |
|---|
| 552 | zum(1)= zu(1) |
|---|
| 553 | zvm(1)= zv(1) |
|---|
| 554 | ztm(1) = ztc(1) |
|---|
| 555 | else |
|---|
| 556 | ! Surface sublimation: |
|---|
| 557 | ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep |
|---|
| 558 | zum(1) = 0 |
|---|
| 559 | zvm(1) = 0 |
|---|
| 560 | end if |
|---|
| 561 | do iq=1,nq |
|---|
| 562 | zqm(1,iq)=0. ! most tracer do not condense ! |
|---|
| 563 | enddo |
|---|
| 564 | ! Special case if the tracer is n2 gas |
|---|
| 565 | if (igcm_n2.ne.0) zqm(1,igcm_n2)=1. |
|---|
| 566 | |
|---|
| 567 | !!! Source haze: 0.02 pourcent when n2 sublimes |
|---|
| 568 | IF (source_haze) THEN |
|---|
| 569 | IF (pdicen2(ig).lt.0) THEN |
|---|
| 570 | DO iq=1,nq |
|---|
| 571 | tname=noms(iq) |
|---|
| 572 | if (tname(1:4).eq."haze") then |
|---|
| 573 | !zqm(1,iq)=0.02 |
|---|
| 574 | !zqm(1,iq)=-pdicen2(ig)*0.02 |
|---|
| 575 | zqm(1,iq)=-pdicen2(ig)*ptimestep*0.02 |
|---|
| 576 | !zqm(10,iq)=-pdicen2(ig)*ptimestep*100. |
|---|
| 577 | !zqm(1,iq)=-pdicen2(ig)*1000000. |
|---|
| 578 | |
|---|
| 579 | endif |
|---|
| 580 | ENDDO |
|---|
| 581 | ENDIF |
|---|
| 582 | ENDIF |
|---|
| 583 | ztm(klev+1)= ztc(klev) ! should not be used, but... |
|---|
| 584 | zum(klev+1)= zu(klev) ! should not be used, but... |
|---|
| 585 | zvm(klev+1)= zv(klev) ! should not be used, but... |
|---|
| 586 | do iq=1,nq |
|---|
| 587 | zqm(klev+1,iq)= zq(klev,iq) |
|---|
| 588 | enddo |
|---|
| 589 | |
|---|
| 590 | ! Tendencies on T, U, V, Q |
|---|
| 591 | ! """"""""""""""""""""""" |
|---|
| 592 | DO l=1,klev |
|---|
| 593 | |
|---|
| 594 | ! Tendencies on T |
|---|
| 595 | zdtsig(ig,l) = (1/masse(l)) * & |
|---|
| 596 | ( zmflux(l)*(ztm(l) - ztc(l)) & |
|---|
| 597 | - zmflux(l+1)*(ztm(l+1) - ztc(l)) & |
|---|
| 598 | + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) ) |
|---|
| 599 | |
|---|
| 600 | ! Tendencies on U |
|---|
| 601 | pduc(ig,l) = (1/masse(l)) * & |
|---|
| 602 | ( zmflux(l)*(zum(l) - zu(l))& |
|---|
| 603 | - zmflux(l+1)*(zum(l+1) - zu(l)) ) |
|---|
| 604 | |
|---|
| 605 | ! Tendencies on V |
|---|
| 606 | pdvc(ig,l) = (1/masse(l)) * & |
|---|
| 607 | ( zmflux(l)*(zvm(l) - zv(l)) & |
|---|
| 608 | - zmflux(l+1)*(zvm(l+1) - zv(l)) ) |
|---|
| 609 | |
|---|
| 610 | END DO |
|---|
| 611 | |
|---|
| 612 | ! Tendencies on Q |
|---|
| 613 | do iq=1,nq |
|---|
| 614 | if (iq.eq.igcm_n2) then |
|---|
| 615 | ! SPECIAL Case when the tracer IS N2 : |
|---|
| 616 | DO l=1,klev |
|---|
| 617 | pdqc(ig,l,iq)= (1/masse(l)) * & |
|---|
| 618 | ( zmflux(l)*(zqm(l,iq) - zq(l,iq)) & |
|---|
| 619 | - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))& |
|---|
| 620 | + zcondicea(ig,l)*(zq(l,iq)-1.) ) |
|---|
| 621 | END DO |
|---|
| 622 | else |
|---|
| 623 | DO l=1,klev |
|---|
| 624 | pdqc(ig,l,iq)= (1/masse(l)) * & |
|---|
| 625 | ( zmflux(l)*(zqm(l,iq) - zq(l,iq)) & |
|---|
| 626 | - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) & |
|---|
| 627 | + zcondicea(ig,l)*zq(l,iq) ) |
|---|
| 628 | END DO |
|---|
| 629 | end if |
|---|
| 630 | enddo |
|---|
| 631 | ! Update variables at the end of each subtimestep. |
|---|
| 632 | do l=1,klev |
|---|
| 633 | ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep |
|---|
| 634 | zu(l) =zu(l) + pduc(ig,l) *subtimestep |
|---|
| 635 | zv(l) =zv(l) + pdvc(ig,l) *subtimestep |
|---|
| 636 | do iq=1,nq |
|---|
| 637 | zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep |
|---|
| 638 | enddo |
|---|
| 639 | end do |
|---|
| 640 | enddo ! loop on nsubtimestep |
|---|
| 641 | ! Recomputing Total tendencies |
|---|
| 642 | do l=1,klev |
|---|
| 643 | pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep |
|---|
| 644 | pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep |
|---|
| 645 | pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep |
|---|
| 646 | do iq=1,nq |
|---|
| 647 | pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep |
|---|
| 648 | |
|---|
| 649 | |
|---|
| 650 | ! Correction temporaire |
|---|
| 651 | if (iq.eq.igcm_n2) then |
|---|
| 652 | if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) & |
|---|
| 653 | .lt.0.01) then ! if n2 < 1 % ! |
|---|
| 654 | write(*,*) 'Warning: n2 < 1%' |
|---|
| 655 | pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq) |
|---|
| 656 | end if |
|---|
| 657 | end if |
|---|
| 658 | |
|---|
| 659 | enddo |
|---|
| 660 | end do |
|---|
| 661 | ! *******************************TEMPORAIRE ****************** |
|---|
| 662 | if (klon.eq.1) then |
|---|
| 663 | write(*,*) 'nsubtimestep=' ,nsubtimestep |
|---|
| 664 | write(*,*) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g |
|---|
| 665 | write(*,*) 'masse apres' , masse(1) |
|---|
| 666 | write(*,*) 'zmflux*DT, l=1' , zmflux(1)*ptimestep |
|---|
| 667 | write(*,*) 'zmflux*DT, l=2' , zmflux(2)*ptimestep |
|---|
| 668 | write(*,*) 'pq, l=1,2,3' , pq(1,1,1), pq(1,2,1),pq(1,3,1) |
|---|
| 669 | write(*,*) 'zq, l=1,2,3' , zq(1,1), zq(2,1),zq(3,1) |
|---|
| 670 | write(*,*) 'dq*Dt l=1' , pdq(1,1,1)*ptimestep |
|---|
| 671 | write(*,*) 'dqcond*Dt l=1' , pdqc(1,1,1)*ptimestep |
|---|
| 672 | end if |
|---|
| 673 | |
|---|
| 674 | ! *********************************************************** |
|---|
| 675 | end if ! if (condsub) |
|---|
| 676 | END DO ! loop on ig |
|---|
| 677 | endif ! not fast |
|---|
| 678 | |
|---|
| 679 | ! ************ Option Olkin to prevent N2 effect in the south******** |
|---|
| 680 | 112 continue |
|---|
| 681 | if (olkin) then |
|---|
| 682 | DO ig=1,klon |
|---|
| 683 | if (latitude(ig).lt.0) then |
|---|
| 684 | pdtsrfc(ig) = max(0.,pdtsrfc(ig)) |
|---|
| 685 | pdpsrf(ig) = 0. |
|---|
| 686 | pdicen2(ig) = 0. |
|---|
| 687 | do l=1,klev |
|---|
| 688 | pdtc(ig,l) = max(0.,zdtlatent(ig,l)) |
|---|
| 689 | pduc(ig,l) = 0. |
|---|
| 690 | pdvc(ig,l) = 0. |
|---|
| 691 | do iq=1,nq |
|---|
| 692 | pdqc(ig,l,iq) = 0 |
|---|
| 693 | enddo |
|---|
| 694 | end do |
|---|
| 695 | end if |
|---|
| 696 | END DO |
|---|
| 697 | end if |
|---|
| 698 | ! ******************************************************************* |
|---|
| 699 | |
|---|
| 700 | ! *************************************************************** |
|---|
| 701 | ! Ecriture des diagnostiques |
|---|
| 702 | ! *************************************************************** |
|---|
| 703 | |
|---|
| 704 | ! DO l=1,klev |
|---|
| 705 | ! DO ig=1,klon |
|---|
| 706 | ! Taux de cond en kg.m-2.pa-1.s-1 |
|---|
| 707 | ! tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1)) |
|---|
| 708 | ! Taux de cond en kg.m-3.s-1 |
|---|
| 709 | ! tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l)) |
|---|
| 710 | ! END DO |
|---|
| 711 | ! END DO |
|---|
| 712 | ! call WRITEDIAGFI(klon,'tconda1', & |
|---|
| 713 | ! 'Taux de condensation N2 atmospherique /Pa', & |
|---|
| 714 | ! 'kg.m-2.Pa-1.s-1',3,tconda1) |
|---|
| 715 | ! call WRITEDIAGFI(klon,'tconda2', & |
|---|
| 716 | ! 'Taux de condensation N2 atmospherique /m', & |
|---|
| 717 | ! 'kg.m-3.s-1',3,tconda2) |
|---|
| 718 | |
|---|
| 719 | |
|---|
| 720 | return |
|---|
| 721 | end subroutine condense_n2 |
|---|
| 722 | |
|---|
| 723 | !------------------------------------------------------------------------- |
|---|
| 724 | |
|---|
| 725 | real function tcond_n2(p,vmr) |
|---|
| 726 | ! Calculates the condensation temperature for N2 at pressure P and vmr |
|---|
| 727 | implicit none |
|---|
| 728 | real, intent(in):: p,vmr |
|---|
| 729 | |
|---|
| 730 | ! tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr)) |
|---|
| 731 | IF (p.ge.0.529995) then |
|---|
| 732 | ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB |
|---|
| 733 | ! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e5*0.98)*log(1./(0.125570*1.e5)*p*vmr)) |
|---|
| 734 | tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr)) |
|---|
| 735 | ELSE |
|---|
| 736 | ! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT |
|---|
| 737 | ! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e5*1.09)*log(1./(0.508059)*p*vmr)) |
|---|
| 738 | tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr)) |
|---|
| 739 | ENDIF |
|---|
| 740 | return |
|---|
| 741 | end function tcond_n2 |
|---|
| 742 | |
|---|
| 743 | !------------------------------------------------------------------------- |
|---|
| 744 | |
|---|
| 745 | real function pcond_n2(t,vmr) |
|---|
| 746 | ! Calculates the condensation pressure for N2 at temperature T and vmr |
|---|
| 747 | implicit none |
|---|
| 748 | real, intent(in):: t,vmr |
|---|
| 749 | |
|---|
| 750 | ! tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr)) |
|---|
| 751 | IF (t.ge.35.6) then |
|---|
| 752 | ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB |
|---|
| 753 | ! pcond_n2 = 0.125570*1.e5/vmr*exp((2.5e5*0.98)/296.925*(1./63.1470-1./t)) |
|---|
| 754 | pcond_n2 = 0.125570e5/vmr*exp(825.1241896*(1./63.147-1./t)) |
|---|
| 755 | ELSE |
|---|
| 756 | ! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB |
|---|
| 757 | ! pcond_n2 = 0.508059/vmr*exp((2.5e5*1.09)/296.925*(1./35.6-1./t)) |
|---|
| 758 | pcond_n2 = 0.508059/vmr*exp(917.7401701*(1./35.6-1./t)) |
|---|
| 759 | ENDIF |
|---|
| 760 | return |
|---|
| 761 | end function pcond_n2 |
|---|
| 762 | |
|---|
| 763 | !------------------------------------------------------------------------- |
|---|
| 764 | |
|---|
| 765 | real function glob_average2d(var) |
|---|
| 766 | ! Calculates the global average of variable var |
|---|
| 767 | use comgeomfi_h |
|---|
| 768 | use dimphy, only: klon |
|---|
| 769 | USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat |
|---|
| 770 | use geometry_mod, only: cell_area, latitude |
|---|
| 771 | |
|---|
| 772 | implicit none |
|---|
| 773 | |
|---|
| 774 | ! INTEGER klon |
|---|
| 775 | REAL var(klon) |
|---|
| 776 | INTEGER ig |
|---|
| 777 | |
|---|
| 778 | glob_average2d = 0. |
|---|
| 779 | DO ig=2,klon-1 |
|---|
| 780 | glob_average2d = glob_average2d + var(ig)*cell_area(ig) |
|---|
| 781 | END DO |
|---|
| 782 | glob_average2d = glob_average2d + var(1)*cell_area(1)*nbp_lon |
|---|
| 783 | glob_average2d = glob_average2d + var(klon)*cell_area(klon)*nbp_lon |
|---|
| 784 | glob_average2d = glob_average2d/(totarea+(cell_area(1)+cell_area(klon))*(nbp_lon-1)) |
|---|
| 785 | |
|---|
| 786 | end function glob_average2d |
|---|
| 787 | |
|---|
| 788 | ! ***************************************************************** |
|---|
| 789 | |
|---|
| 790 | subroutine vl1d(klev,q,pente_max,masse,w,qm) |
|---|
| 791 | ! |
|---|
| 792 | ! Operateur de moyenne inter-couche pour calcul de transport type |
|---|
| 793 | ! Van-Leer " pseudo amont " dans la verticale |
|---|
| 794 | ! q,w sont des arguments d'entree pour le s-pg .... |
|---|
| 795 | ! masse : masse de la couche Dp/g |
|---|
| 796 | ! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2) |
|---|
| 797 | ! pente_max = 2 conseillee |
|---|
| 798 | ! -------------------------------------------------------------------- |
|---|
| 799 | IMPLICIT NONE |
|---|
| 800 | |
|---|
| 801 | ! Arguments: |
|---|
| 802 | ! ---------- |
|---|
| 803 | integer klev |
|---|
| 804 | real masse(klev),pente_max |
|---|
| 805 | REAL q(klev),qm(klev+1) |
|---|
| 806 | REAL w(klev+1) |
|---|
| 807 | ! |
|---|
| 808 | ! Local |
|---|
| 809 | ! --------- |
|---|
| 810 | ! |
|---|
| 811 | INTEGER l |
|---|
| 812 | ! |
|---|
| 813 | real dzq(klev),dzqw(klev),adzqw(klev),dzqmax |
|---|
| 814 | real sigw, Mtot, MQtot |
|---|
| 815 | integer m |
|---|
| 816 | |
|---|
| 817 | |
|---|
| 818 | ! On oriente tout dans le sens de la pression |
|---|
| 819 | ! W > 0 WHEN DOWN !!!!!!!!!!!!! |
|---|
| 820 | |
|---|
| 821 | do l=2,klev |
|---|
| 822 | dzqw(l)=q(l-1)-q(l) |
|---|
| 823 | adzqw(l)=abs(dzqw(l)) |
|---|
| 824 | enddo |
|---|
| 825 | |
|---|
| 826 | do l=2,klev-1 |
|---|
| 827 | if(dzqw(l)*dzqw(l+1).gt.0.) then |
|---|
| 828 | dzq(l)=0.5*(dzqw(l)+dzqw(l+1)) |
|---|
| 829 | else |
|---|
| 830 | dzq(l)=0. |
|---|
| 831 | endif |
|---|
| 832 | dzqmax=pente_max*min(adzqw(l),adzqw(l+1)) |
|---|
| 833 | dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l)) |
|---|
| 834 | enddo |
|---|
| 835 | |
|---|
| 836 | dzq(1)=0. |
|---|
| 837 | dzq(klev)=0. |
|---|
| 838 | |
|---|
| 839 | do l = 1,klev-1 |
|---|
| 840 | |
|---|
| 841 | ! Regular scheme (transfered mass < layer mass) |
|---|
| 842 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 843 | if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then |
|---|
| 844 | sigw=w(l+1)/masse(l+1) |
|---|
| 845 | qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1)) |
|---|
| 846 | else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then |
|---|
| 847 | sigw=w(l+1)/masse(l) |
|---|
| 848 | qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l)) |
|---|
| 849 | |
|---|
| 850 | ! Extended scheme (transfered mass > layer mass) |
|---|
| 851 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 852 | else if(w(l+1).gt.0.) then |
|---|
| 853 | m=l+1 |
|---|
| 854 | Mtot = masse(m) |
|---|
| 855 | MQtot = masse(m)*q(m) |
|---|
| 856 | !if (m.lt.klev) then ! because some compilers will have problems |
|---|
| 857 | ! ! evaluating masse(klev+1) |
|---|
| 858 | do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1)))) |
|---|
| 859 | m=m+1 |
|---|
| 860 | Mtot = Mtot + masse(m) |
|---|
| 861 | MQtot = MQtot + masse(m)*q(m) |
|---|
| 862 | ! if (m.eq.klev) exit |
|---|
| 863 | end do |
|---|
| 864 | !endif |
|---|
| 865 | if (m.lt.klev) then |
|---|
| 866 | sigw=(w(l+1)-Mtot)/masse(m+1) |
|---|
| 867 | qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)* & |
|---|
| 868 | (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) ) |
|---|
| 869 | else |
|---|
| 870 | w(l+1) = Mtot |
|---|
| 871 | qm(l+1) = Mqtot / Mtot |
|---|
| 872 | write(*,*) 'top layer is disapearing !' |
|---|
| 873 | stop |
|---|
| 874 | end if |
|---|
| 875 | else ! if(w(l+1).lt.0) |
|---|
| 876 | m = l-1 |
|---|
| 877 | Mtot = masse(m+1) |
|---|
| 878 | MQtot = masse(m+1)*q(m+1) |
|---|
| 879 | if (m.gt.0) then ! because some compilers will have problems |
|---|
| 880 | ! evaluating masse(0) |
|---|
| 881 | do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m)))) |
|---|
| 882 | m=m-1 |
|---|
| 883 | Mtot = Mtot + masse(m+1) |
|---|
| 884 | MQtot = MQtot + masse(m+1)*q(m+1) |
|---|
| 885 | if (m.eq.0) exit |
|---|
| 886 | end do |
|---|
| 887 | endif |
|---|
| 888 | if (m.gt.0) then |
|---|
| 889 | sigw=(w(l+1)+Mtot)/masse(m) |
|---|
| 890 | qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* & |
|---|
| 891 | (q(m)-0.5*(1.+sigw)*dzq(m)) ) |
|---|
| 892 | else |
|---|
| 893 | qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1)) |
|---|
| 894 | end if |
|---|
| 895 | end if |
|---|
| 896 | enddo |
|---|
| 897 | |
|---|
| 898 | return |
|---|
| 899 | end subroutine vl1d |
|---|
| 900 | |
|---|
