Changeset 665
- Timestamp:
- May 23, 2012, 10:33:16 AM (13 years ago)
- Location:
- trunk
- Files:
-
- 9 edited
-
MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90 (modified) (2 diffs)
-
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM (modified) (1 diff)
-
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/aeropacity_use.F (modified) (4 diffs)
-
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/callradite_use.F (modified) (2 diffs)
-
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/initracer_use.F (modified) (8 diffs)
-
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/physiq_modif_use.F (modified) (59 diffs)
-
MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F (modified) (2 diffs)
-
MESOSCALE_DEV/MANUAL/SRC/faq.tex (modified) (1 diff)
-
UTIL/PYTHON/planetoplot.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90
r549 r665 662 662 if (ierr.ne.NF_NOERR) then 663 663 write(*,*) "...No ccn - CCN mass set to 0" 664 dustfile(:,:,:,:)=0.664 ccnfile(:,:,:,:)=0. 665 665 else 666 666 ierr=NF_GET_VAR_REAL(nid,nvarid,ccnfile) … … 671 671 if (ierr.ne.NF_NOERR) then 672 672 write(*,*) "...No ccnN - CCN number set to 0" 673 dustnfile(:,:,:,:)=0.673 ccnnfile(:,:,:,:)=0. 674 674 else 675 675 ierr=NF_GET_VAR_REAL(nid,nvarid,ccnnfile) -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
r549 r665 87 87 state real THETA ij misc 1 - rd "THETA" "SLOPE INCLINATION" "deg" #SAVEMARS2 theta_sl 88 88 state real PSI ij misc 1 - rd "PSI" "SLOPE ORIENTATION" "deg" #SAVEMARS2 psi_sl 89 state real TAU_DUST ij misc 1 - r d "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" "" #SAVEMARS2 tauref89 state real TAU_DUST ij misc 1 - rhd "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" "" #SAVEMARS2 tauref 90 90 state real SWDOWNZ ij misc 1 - rd "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2" #SAVEMARS2 fluxsurf_sw_tot 91 91 state real LWDOWNZ ij misc 1 - rd "LWDOWNZ" "DOWNWARD LW FLUX AT SURFACE" "W m-2" #SAVEMARS2 fluxsurf_lw 92 92 state real SWUP ij misc 1 - rd "SWUP" "UPWARD SW FLUX AT TOP" "W m-2" #SAVEMARS2 fluxtop_sw_tot 93 93 state real LWUP ij misc 1 - rd "LWUP" "UPWARD LW FLUX AT TOP" "W m-2" #SAVEMARS2 fluxtop_lw 94 state real MTOT ij misc 1 - r d "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic" #SAVEMARS2 mtot95 state real ICETOT ij misc 1 - r d "ICETOT" "TOTAL MASS WATER ICE" "kg m-2" #SAVEMARS2 icetot96 state real RAVE ij misc 1 - r d "RAVE" "MEAN ICE RADIUS" "m" #SAVEMARS2 rave94 state real MTOT ij misc 1 - rhd "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic" #SAVEMARS2 mtot 95 state real ICETOT ij misc 1 - rhd "ICETOT" "TOTAL MASS WATER ICE" "kg m-2" #SAVEMARS2 icetot 96 state real RAVE ij misc 1 - rhd "RAVE" "MEAN ICE RADIUS" "m" #SAVEMARS2 rave 97 97 state real RICE ikj misc 1 - rd "RICE" "ICE RADIUS" "m" #SAVEMARS3 rice 98 98 state real HR_SW ikj misc 1 - rd "HR_SW" "HEATING RATE SW" "K/s" #SAVEMARS3 zdtsw 99 99 state real HR_LW ikj misc 1 - rd "HR_LW" "HEATING RATE LW" "K/s" #SAVEMARS3 zdtlw 100 100 state real HR_SH ikj misc 1 - rd "HR_SH" "HEATING RATE sens. heat" "K/s" #SAVEMARS3 zdtdif 101 state real QSURFICE ij misc 1 - r d "QSURFICE" "WATER ICE AT SURFACE" "kg m-2" #SAVEMARS2 qsurfice102 state real RDUST ikj misc 1 - r d "RDUST" "DUST RADIUS" "m" #SAVEMARS3 rdust101 state real QSURFICE ij misc 1 - rhd "QSURFICE" "WATER ICE AT SURFACE" "kg m-2" #SAVEMARS2 qsurfice 102 state real RDUST ikj misc 1 - rhd "RDUST" "DUST RADIUS" "m" #SAVEMARS3 rdust 103 103 state real HR_NIR ikj misc 1 - rd "HR_NIR" "HEATING RATE nirco2" "K/s" #SAVEMARS3 zdtnirco2 104 104 state real HR_NLTE ikj misc 1 - rd "HR_NLTE" "HEATING RATE nlte" "K/s" #SAVEMARS3 zdtnlte 105 state real ALBBARE ij misc 1 - r d "ALBBARE" "SOIL ALBEDO" "" #SAVEMARS2 albedodat106 state real VMR_ICE ikj misc 1 - r d "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm" #SAVEMARS3 vmr107 state real TAU_ICE ij misc 1 - r d "TAU_ICE" "CLOUD OD at 825 cm-1 TES" "" #SAVEMARS2 tauTES105 state real ALBBARE ij misc 1 - rhd "ALBBARE" "SOIL ALBEDO" "" #SAVEMARS2 albedodat 106 state real VMR_ICE ikj misc 1 - rhd "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm" #SAVEMARS3 vmr 107 state real TAU_ICE ij misc 1 - rhd "TAU_ICE" "CLOUD OD at 825 cm-1 TES" "" #SAVEMARS2 tauTES 108 108 state real PDTZ ikj misc 1 - rd "PDT" "TEMP TENDENCY" "K s-1" #SAVEMARS3 pdt 109 state real ZMAX_TH ij misc 1 - r d "ZMAX_TH" "MAXIMUM LEVEL REACHED IN TH" "m" #SAVEMARS2 zmax_th109 state real ZMAX_TH ij misc 1 - rhd "ZMAX_TH" "MAXIMUM LEVEL REACHED IN TH" "m" #SAVEMARS2 zmax_th 110 110 state real HFMAX_TH ij misc 1 - rd "HFMAX_TH" "MAXIMUM TH HEAT FLUX" "m.K/s" #SAVEMARS2 hfmax_th 111 state real WSTAR ij misc 1 - r d "WSTAR" "FREE CONVECTION VELOCITY FROM TH" "m/s" #SAVEMARS2 wstar111 state real WSTAR ij misc 1 - rhd "WSTAR" "FREE CONVECTION VELOCITY FROM TH" "m/s" #SAVEMARS2 wstar 112 112 state real Z0SET ij misc 1 - rd "Z0SET" "SET SURFACE ROUGHNESS" "m" #SAVEMARS2 z0 113 113 -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/aeropacity_use.F
r308 r665 189 189 msolsir(1:nlayer,1:naerkind)=0 190 190 mqextsqabs(1:nlayer,1:naerkind)=0 191 WRITE(*,*) "Typical profiles of solsir and Qext/Qabs(IR):" 191 WRITE(*,*) "Typical profiles of Qext(vis)/Qext(IR)" 192 WRITE(*,*) " and Qext(IR)/Qabs(IR):" 192 193 DO iaer = 1, naerkind ! Loop on aerosol kind 193 194 WRITE(*,*) "Aerosol # ",iaer … … 210 211 ! otherwise default value read from starfi.nc file will be used) 211 212 call getin("tauvis",tauvis) 213 214 ! Some information about the water cycle 215 IF (water) THEN 216 write(*,*) "water_param CCN reduc. fac. ", ccn_factor 217 ENDIF 212 218 213 219 firstcall=.false. … … 445 451 taudusttmp(ig) = taudusttmp(ig) + 446 452 & aerosol(ig,l,iaerdust(iaer)) 447 448 453 ENDDO 449 454 ENDDO … … 590 595 ! ccn(ig,l) = ccn(ig,l) / 8. 591 596 c ----------------------------------------------------------------- 592 write(*,*) "water_param CCN reduc. fac. ", ccn_factor593 597 DO l=1,nlayer 594 598 DO ig=1,ngrid -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/callradite_use.F
r308 r665 264 264 c PLEASE MAKE SURE that you set up the right number of 265 265 c scatterers in dimradmars.h (naerkind); 266 c name_iaer(1) = "dust_conrath" !! poussiere classique 267 name_iaer(1) = "dust_doubleq" 268 cc name_iaer(2) = "dust_submicron" !! JB: experimental 269 c name_iaer(2) = "h2o_ice" 266 name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile 267 IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme 268 if (nqmx.gt.1) then 269 ! trick to avoid problems compiling with 1 tracer 270 ! and picky compilers who know name_iaer(2) is out of bounds 271 j=2 272 IF (water.AND.activice) name_iaer(j) = "h2o_ice" !! radiatively-active clouds 273 IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff 274 endif ! of if (nqmx.gt.1) 270 275 c ---------------------------------------------------------- 271 276 … … 389 394 c Computing aerosol optical depth in each layer: 390 395 CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls, 391 & pq,ccn,tauref,tau,aerosol,reffrad,nueffrad, QREFvis3d, 392 & QREFir3d,omegaREFvis3d,omegaREFir3d,zdqnorm) 396 & pq,ccn,tauref,tau,aerosol,reffrad,nueffrad, QREFvis3d, 397 & QREFir3d,omegaREFvis3d,omegaREFir3d,zdqnorm) 398 393 399 c Starting loop on sub-domain 394 400 c ---------------------------- -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/initracer_use.F
r308 r665 43 43 #include "advtrac.h" 44 44 #include "comgeomfi.h" 45 #include "watercap.h"46 45 #include "chimiedata.h" 47 46 … … 111 110 write(*,*) 'initracer: error expected dustbin=2' 112 111 else 113 noms(1)='dust_mass' ! dust mass mixing ratio 114 noms(2)='dust_number' ! dust number mixing ratio 112 ! noms(1)='dust_mass' ! dust mass mixing ratio 113 ! noms(2)='dust_number' ! dust number mixing ratio 114 ! same as above, but avoid explict possible out-of-bounds on array 115 noms(1)='dust_mass' ! dust mass mixing ratio 116 do iq=2,2 117 noms(iq)='dust_number' ! dust number mixing ratio 118 enddo 115 119 endif 116 120 endif … … 120 124 noms(nqmx)='h2o_vap' 121 125 mmol(nqmx)=18. 122 noms(nqmx-1)='h2o_ice' 123 mmol(nqmx-1)=18. 126 ! noms(nqmx-1)='h2o_ice' 127 ! mmol(nqmx-1)=18. 128 !"loop" to avoid potential out-of-bounds in array 129 do iq=nqmx-1,nqmx-1 130 noms(iq)='h2o_ice' 131 mmol(iq)=18. 132 enddo 124 133 endif 125 134 ! 3. Chemistry … … 157 166 if (oldnames.and.water) then 158 167 write(*,*)'initracer: moving surface water ice to index ',nqmx-1 159 qsurf(1:ngridmx,nqmx-1)=qsurf(1:ngridmx,nqmx) 168 ! "loop" to avoid potential out-of-bounds on arrays 169 do iq=nqmx-1,nqmx-1 170 qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1) 171 enddo 160 172 qsurf(1:ngridmx,nqmx)=0 161 173 endif … … 218 230 igcm_ar=0 219 231 igcm_ar_n2=0 232 igcm_n=0 233 igcm_no=0 234 igcm_no2=0 235 igcm_n2d=0 236 igcm_co2plus=0 237 igcm_oplus=0 238 igcm_o2plus=0 239 igcm_coplus=0 240 igcm_cplus=0 241 igcm_nplus=0 242 igcm_noplus=0 243 igcm_n2plus=0 244 igcm_hplus=0 245 igcm_elec=0 246 220 247 221 248 ! 1. find dust tracers … … 320 347 count=count+1 321 348 endif 349 if (noms(iq).eq."n") then 350 igcm_n=iq 351 mmol(igcm_n)=14. 352 count=count+1 353 endif 354 if (noms(iq).eq."no") then 355 igcm_no=iq 356 mmol(igcm_no)=30. 357 count=count+1 358 endif 359 if (noms(iq).eq."no2") then 360 igcm_no2=iq 361 mmol(igcm_no2)=46. 362 count=count+1 363 endif 364 if (noms(iq).eq."n2d") then 365 igcm_n2d=iq 366 mmol(igcm_n2d)=28. 367 count=count+1 368 endif 369 if (noms(iq).eq."co2plus") then 370 igcm_co2plus=iq 371 mmol(igcm_co2plus)=44. 372 count=count+1 373 endif 374 if (noms(iq).eq."oplus") then 375 igcm_oplus=iq 376 mmol(igcm_oplus)=16. 377 count=count+1 378 endif 379 if (noms(iq).eq."o2plus") then 380 igcm_o2plus=iq 381 mmol(igcm_o2plus)=32. 382 count=count+1 383 endif 384 if (noms(iq).eq."coplus") then 385 igcm_coplus=iq 386 mmol(igcm_coplus)=28. 387 count=count+1 388 endif 389 if (noms(iq).eq."cplus") then 390 igcm_cplus=iq 391 mmol(igcm_cplus)=12. 392 count=count+1 393 endif 394 if (noms(iq).eq."nplus") then 395 igcm_nplus=iq 396 mmol(igcm_nplus)=14. 397 count=count+1 398 endif 399 if (noms(iq).eq."noplus") then 400 igcm_noplus=iq 401 mmol(igcm_noplus)=30. 402 count=count+1 403 endif 404 if (noms(iq).eq."n2plus") then 405 igcm_n2plus=iq 406 mmol(igcm_n2plus)=28. 407 count=count+1 408 endif 409 if (noms(iq).eq."hplus") then 410 igcm_hplus=iq 411 mmol(igcm_hplus)=1. 412 count=count+1 413 endif 414 if (noms(iq).eq."elec") then 415 igcm_elec=iq 416 mmol(igcm_elec)=1./1822.89 417 count=count+1 418 endif 322 419 if (noms(iq).eq."h2o_vap") then 323 420 igcm_h2o_vap=iq … … 336 433 count=count+1 337 434 endif 435 436 338 437 enddo ! of do iq=1,nqmx 339 438 ! count=count+nbqchem … … 486 585 alpha_lift(igcm_h2o_vap) =0. 487 586 alpha_devil(igcm_h2o_vap)=0. 488 489 c "Dryness coefficient" controlling the evaporation and490 c sublimation from the ground water ice (close to 1)491 c HERE, the goal is to correct for the fact492 c that the simulated permanent water ice polar caps493 c is larger than the actual cap and the atmospheric494 c opacity not always realistic.495 496 do ig=1,ngridmx497 if (ngridmx.ne.1) watercaptag(ig)=.false.498 dryness(ig) = 1.499 enddo500 501 IF (caps) THEN502 c Perennial H20 north cap defined by watercaptag=true (allows surface to be503 c hollowed by sublimation in vdifc).504 yeyey = 0505 do ig=1,ngridmx506 ! !!! TESTS TESTS outliers507 ! !!! TESTS TESTS outliers508 ! if ( ( lati(ig)*180./pi .ge. 75 ) .and.509 ! . ( lati(ig)*180./pi .le. 77 ) .and.510 ! . ( ( ( long(ig)*180./pi .ge. 000. ) .and.511 ! . ( long(ig)*180./pi .le. 120. ) )512 ! . .or.513 ! . ( ( long(ig)*180./pi .ge. -130. ) .and.514 ! . ( long(ig)*180./pi .le. -115. ) ) ) ) then515 ! if (yeyey .eq. 0) then !!! 1/2 en 64x48 sinon trop large en lat516 ! write(*,*) "outliers ", lati(ig)*180./pi, long(ig)*180./pi517 ! if (ngridmx.ne.1) watercaptag(ig)=.true.518 ! dryness(ig) = 1.519 ! albedodat(ig) = 0.45 !! comme alb_surfice520 ! yeyey = 1521 ! else522 ! yeyey = 0523 ! endif524 ! endif525 ! !!! TESTS TESTS outliers526 ! !!! TESTS TESTS outliers527 !528 ! !!! TESTS TESTS addcap529 ! !!! TESTS TESTS addcap530 ! if ( ( lati(ig)*180./pi .ge. 82 ) .and.531 ! . ( lati(ig)*180./pi .le. 84 ) .and.532 ! . ( ( long(ig)*180./pi .gt. -030. ) .and.533 ! . ( long(ig)*180./pi .lt. 090. ) ) ) then534 ! write(*,*) "capadd ", lati(ig)*180./pi, long(ig)*180./pi535 ! if (ngridmx.ne.1) watercaptag(ig)=.true.536 ! albedodat(ig) = 0.45 !! comme alb_surfice537 ! dryness(ig) = 1.538 ! endif539 ! !!! TESTS TESTS addcap540 ! !!! TESTS TESTS addcap541 542 if (lati(ig)*180./pi.gt.84) then543 if (ngridmx.ne.1) watercaptag(ig)=.true.544 dryness(ig) = 1.545 c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)546 c if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then547 c if (ngridmx.ne.1) watercaptag(ig)=.true.548 c dryness(ig) = 1.549 c endif550 c if (lati(ig)*180./pi.ge.85) then551 c if (ngridmx.ne.1) watercaptag(ig)=.true.552 c dryness(ig) = 1.553 c endif554 endif ! (lati>80 deg)555 end do ! (ngridmx)556 ENDIF ! (caps)557 558 587 if(water.and.(nqmx.ge.2)) then 559 588 radius(igcm_h2o_ice)=3.e-6 -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/physiq_modif_use.F
r308 r665 1 SUBROUTINE physiq(ngrid,nlayer,nq, 2 $ firstcall,lastcall, 3 $ pday,ptime,ptimestep, 4 $ pplev,pplay,pphi, 5 $ pu,pv,pt,pq, 6 $ pw, 7 $ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn) 8 1 SUBROUTINE physiq( 2 $ ngrid,nlayer,nq 3 $ ,firstcall,lastcall 4 $ ,pday,ptime,ptimestep 5 $ ,pplev,pplay,pphi 6 $ ,pu,pv,pt,pq 7 $ ,pw 8 $ ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn 9 #ifdef MESOSCALE 10 #include "meso_inc/meso_inc_invar.F" 11 #endif 12 $ ) 9 13 10 14 IMPLICIT NONE … … 42 46 c 8. Contribution to tendencies due to thermosphere 43 47 c 9. Surface and sub-surface temperature calculations 44 c 10. Renormalisation 45 c 11. Write outputs : 48 c 10. Write outputs : 46 49 c - "startfi", "histfi" (if it's time) 47 50 c - Saving statistics (if "callstats = .true.") 48 51 c - Dumping eof (if "calleofdump = .true.") 49 52 c - Output any needed variables in "diagfi" 50 c 1 2. Diagnostic: mass conservation of tracers53 c 11. Diagnostic: mass conservation of tracers 51 54 c 52 55 c author: … … 61 64 c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009) 62 65 c Nb: See callradite.F for more information. 66 c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags 63 67 c 64 68 c arguments: … … 124 128 125 129 #include "chimiedata.h" 126 #include "watercap.h"127 130 #include "param.h" 128 131 #include "param_v3.h" … … 131 134 #include "netcdf.inc" 132 135 133 136 #include "slope.h" 137 138 #ifdef MESOSCALE 139 #include "wrf_output_2d.h" 140 #include "wrf_output_3d.h" 141 #include "advtrac.h" !!! this is necessary for tracers (in dyn3d) 142 #include "meso_inc/meso_inc_var.F" 143 #endif 134 144 135 145 c Arguments : … … 180 190 REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2) 181 191 REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2) 182 REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy 183 INTEGER ig_vl1 ! Grid Point near VL1 (for diagnostic) 192 REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy 193 194 REAL watercapflag(ngridmx) ! water cap flag 184 195 185 196 c Variables used by the water ice microphysical scheme: … … 187 198 REAL nuice(ngridmx,nlayermx) ! Estimated effective variance 188 199 ! of the size distribution 189 c Albedo of deposited surface ice 190 !!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45 191 REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB 200 201 c Variables used by the slope model 202 REAL sl_ls, sl_lct, sl_lat 203 REAL sl_tau, sl_alb, sl_the, sl_psi 204 REAL sl_fl0, sl_flu 205 REAL sl_ra, sl_di0 206 REAL sky 192 207 193 208 SAVE day_ini, icount … … 195 210 SAVE co2ice,albedo,emis, q2 196 211 SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf 197 SAVE ig_vl1198 212 199 213 REAL stephan … … 251 265 REAL zdqchim(ngridmx,nlayermx,nqmx) 252 266 REAL zdqschim(ngridmx,nqmx) 253 REAL zdqnorm(ngridmx,nlayermx,2) !quantity of dust which have to be added by the dynamical core267 REAL zdqnorm(ngridmx,nlayermx,2) !quantity of dust which have to be added 254 268 255 269 … … 304 318 REAL time_phys 305 319 320 c Variables for PBL 321 322 REAL lmax_th_out(ngridmx),zmax_th(ngridmx) 323 REAL, SAVE :: wmax_th(ngridmx) 324 REAL hfmax_th(ngridmx) 325 REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx) 326 REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx) 327 INTEGER lmax_th(ngridmx) 328 REAL dtke_th(ngridmx,nlayermx+1) 329 REAL zcdv(ngridmx), zcdh(ngridmx) 330 REAL Teta_out(ngridmx),u_out(ngridmx) ! Interpolated teta and u at z_out 331 REAL z_out ! height of interpolation between z0 and z1 332 REAL ustar(ngridmx),tstar(ngridmx) ! friction velocity and friction potential temp 333 REAL zu2(ngridmx) 306 334 c======================================================================= 307 335 … … 315 343 c variables set to 0 316 344 c ~~~~~~~~~~~~~~~~~~ 317 call zerophys(ngrid*nlayer*naerkind,aerosol) 318 call zerophys(ngrid*nlayer,dtrad) 319 call zerophys(ngrid,fluxrad) 345 aerosol(:,:,:)=0 346 dtrad(:,:)=0 347 fluxrad(:)=0 348 349 wmax_th(:)=0. 320 350 321 351 c read startfi 322 352 c ~~~~~~~~~~~~ 323 353 #ifndef MESOSCALE 324 354 ! Read netcdf initial physical parameters. 325 355 CALL phyetat0 ("startfi.nc",0,0, … … 327 357 & day_ini,time_phys, 328 358 & tsurf,tsoil,emis,q2,qsurf,co2ice) 359 #else 360 #include "meso_inc/meso_inc_ini.F" 361 #endif 329 362 330 363 if (pday.ne.day_ini) then … … 337 370 338 371 write (*,*) 'In physiq day_ini =', day_ini 372 373 c initialize tracers 374 c ~~~~~~~~~~~~~~~~~~ 375 tracerdyn=tracer 376 IF (tracer) THEN 377 CALL initracer(qsurf,co2ice) 378 ENDIF ! end tracer 339 379 340 380 c Initialize albedo and orbital calculation … … 358 398 icount=1 359 399 360 361 c initialize tracers 362 c ~~~~~~~~~~~~~~~~~~ 363 tracerdyn=tracer 364 IF (tracer) THEN 365 CALL initracer(qsurf,co2ice) 366 ENDIF ! end tracer 367 368 c Determining gridpoint near Viking Lander 1 (used for diagnostic only) 369 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 370 371 if(ngrid.ne.1) then 372 latvl1= 22.27 373 lonvl1= -47.94 374 ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.) -2 )*iim + 375 & int(1.5+(lonvl1+180)*iim/360.) 376 write(*,*) 'Viking Lander 1 GCM point: lat,lon', 377 & lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi 378 end if 379 400 #ifndef MESOSCALE 380 401 c Initialize thermospheric parameters 381 402 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 382 403 383 404 if (callthermos) call param_read 384 405 #endif 385 406 c Initialize R and Cp as constant 386 407 … … 396 417 397 418 IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN 398 write(*,*)"physiq: water_param Surface ice alb:",alb_surfice 419 write(*,*)"physiq: water_param Surface water ice albedo:", 420 . albedo_h2o_ice 399 421 ENDIF 400 422 … … 415 437 c Initialize various variables 416 438 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 417 call zerophys(ngrid*nlayer, pdv)418 call zerophys(ngrid*nlayer, pdu)419 call zerophys(ngrid*nlayer, pdt)420 call zerophys(ngrid*nlayer*nq, pdq)421 call zerophys(ngrid, pdpsrf)422 call zerophys(ngrid, zflubid)423 call zerophys(ngrid, zdtsurf)424 call zerophys(ngrid*nq, dqsurf)439 pdv(:,:)=0 440 pdu(:,:)=0 441 pdt(:,:)=0 442 pdq(:,:,:)=0 443 pdpsrf(:)=0 444 zflubid(:)=0 445 zdtsurf(:)=0 446 dqsurf(:,:)=0 425 447 igout=ngrid/2+1 426 448 … … 469 491 ENDDO 470 492 493 #ifndef MESOSCALE 471 494 c----------------------------------------------------------------------- 472 495 c 1.2.5 Compute mean mass, cp, and R … … 476 499 call concentrations(pplay,pt,pdt,pq,pdq,ptimestep) 477 500 endif 501 #endif 478 502 c----------------------------------------------------------------------- 479 503 c 2. Compute radiative tendencies : … … 518 542 c Radiative transfer 519 543 c ------------------ 544 520 545 CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, 521 546 $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, … … 523 548 & tauref,tau,aerosol,ccn,rdust,rice,nuice,zdqnorm) 524 549 550 c Outputs for basic check (middle of domain) 551 c ------------------------------------------ 552 print*, 'Ls =',zls*180./pi, 553 & 'check lat lon', lati(igout)*180/pi, 554 & long(igout)*180/pi 555 print*, 'tauref(700 Pa) =',tauref(igout), 556 & ' tau(700 Pa) =',tau(igout,1)*700./pplev(igout,1) 557 558 c --------------------------------------------------------- 559 c Call slope parameterization for direct and scattered flux 560 c --------------------------------------------------------- 561 IF(callslope) THEN 562 print *, 'Slope scheme is on and computing...' 563 DO ig=1,ngrid 564 sl_the = theta_sl(ig) 565 IF (sl_the .ne. 0.) THEN 566 ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2) 567 DO l=1,2 568 sl_lct = ptime*24. + 180.*long(ig)/pi/15. 569 sl_ra = pi*(1.0-sl_lct/12.) 570 sl_lat = 180.*lati(ig)/pi 571 sl_tau = tau(ig,1) 572 sl_alb = albedo(ig,l) 573 sl_psi = psi_sl(ig) 574 sl_fl0 = fluxsurf_sw(ig,l) 575 sl_di0 = 0. 576 if (mu0(ig) .gt. 0.) then 577 sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) 578 sl_di0 = sl_di0*1370./dist_sol/dist_sol 579 sl_di0 = sl_di0/ztim1 580 sl_di0 = fluxsurf_sw(ig,l)*sl_di0 581 endif 582 ! you never know (roundup concern...) 583 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 584 !!!!!!!!!!!!!!!!!!!!!!!!!! 585 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 586 & sl_tau, sl_alb, sl_the, sl_psi, 587 & sl_di0, sl_fl0, sl_flu ) 588 !!!!!!!!!!!!!!!!!!!!!!!!!! 589 fluxsurf_sw(ig,l) = sl_flu 590 ENDDO 591 !!! compute correction on IR flux as well 592 sky= (1.+cos(pi*theta_sl(ig)/180.))/2. 593 fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky 594 ENDIF 595 ENDDO 596 ENDIF 597 525 598 c CO2 near infrared absorption 526 599 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 527 call zerophys(ngrid*nlayer,zdtnirco2)600 zdtnirco2(:,:)=0 528 601 if (callnirco2) then 529 602 call nirco2abs (ngrid,nlayer,pplay,dist_sol, … … 554 627 ENDIF 555 628 556 557 558 629 ENDIF ! of if(mod(icount-1,iradia).eq.0) 559 630 … … 569 640 $ stephan*zplanck(ig)*zplanck(ig) 570 641 fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig) 642 IF(callslope) THEN 643 sky= (1.+cos(pi*theta_sl(ig)/180.))/2. 644 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig) 645 ENDIF 571 646 ENDDO 572 573 647 574 648 DO l=1,nlayer … … 602 676 c 4. Vertical diffusion (turbulent mixing): 603 677 c ----------------------------------------- 604 c 678 605 679 IF (calldifv) THEN 606 607 680 608 681 DO ig=1,ngrid … … 610 683 ENDDO 611 684 612 CALL zerophys(ngrid*nlayer,zdum1)613 CALL zerophys(ngrid*nlayer,zdum2)685 zdum1(:,:)=0 686 zdum2(:,:)=0 614 687 do l=1,nlayer 615 688 do ig=1,ngrid … … 617 690 enddo 618 691 enddo 619 692 693 694 #ifdef MESOSCALE 695 IF (.not.flag_LES) THEN 696 #endif 697 c ---------------------- 698 c Treatment of a special case : using new surface layer (Richardson based) 699 c without using the thermals in gcm and mesoscale can yield problems in 700 c weakly unstable situations when winds are near to 0. For those cases, we add 701 c a unit subgrid gustiness. Remember that thermals should be used we using the 702 c Richardson based surface layer model. 703 IF ( .not.calltherm .and. callrichsl ) THEN 704 DO ig=1, ngridmx 705 IF (zh(ig,1) .lt. tsurf(ig)) THEN 706 wmax_th(ig)=1. 707 ENDIF 708 ENDDO 709 ENDIF 710 c ---------------------- 711 #ifdef MESOSCALE 712 ENDIF 713 #endif 714 715 620 716 c Calling vdif (Martian version WITH CO2 condensation) 621 717 CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk, … … 625 721 $ zdum1,zdum2,zdh,pdq,zflubid, 626 722 $ zdudif,zdvdif,zdhdif,zdtsdif,q2, 627 & zdqdif,zdqsdif) 628 723 & zdqdif,zdqsdif,wmax_th,zcdv,zcdh) 724 725 #ifdef MESOSCALE 726 #include "meso_inc/meso_inc_les.F" 727 #endif 629 728 DO l=1,nlayer 630 729 DO ig=1,ngrid … … 638 737 ENDDO 639 738 640 DO ig=1,ngrid641 zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)642 ENDDO739 DO ig=1,ngrid 740 zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig) 741 ENDDO 643 742 644 743 if (tracer) then … … 662 761 s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) 663 762 ENDDO 763 #ifdef MESOSCALE 764 IF (flag_LES) THEN 765 write(*,*) 'LES mode !' 766 write(*,*) 'Please set calldifv to T in callphys.def' 767 STOP 768 ENDIF 769 #endif 664 770 ENDIF ! of IF (calldifv) 665 771 772 c----------------------------------------------------------------------- 773 c TEST. Thermals : 774 c HIGHLY EXPERIMENTAL, BEWARE !! 775 c ----------------------------- 776 777 if(calltherm) then 778 779 call calltherm_interface(firstcall, 780 $ long,lati,zzlev,zzlay, 781 $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, 782 $ pplay,pplev,pphi,zpopsk, 783 $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th, 784 $ dtke_th,hfmax_th,wmax_th) 785 786 DO l=1,nlayer 787 DO ig=1,ngrid 788 pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l) 789 pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l) 790 pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l) 791 q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep 792 ENDDO 793 ENDDO 794 795 DO ig=1,ngrid 796 q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep 797 ENDDO 798 799 if (tracer) then 800 DO iq=1,nq 801 DO l=1,nlayer 802 DO ig=1,ngrid 803 pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq) 804 ENDDO 805 ENDDO 806 ENDDO 807 endif 808 809 lmax_th_out(:)=real(lmax_th(:)) 810 811 else !of if calltherm 812 lmax_th(:)=0 813 wmax_th(:)=0. 814 lmax_th_out(:)=0. 815 end if 666 816 667 817 c----------------------------------------------------------------------- … … 676 826 ENDDO 677 827 ENDDO 678 CALL zerophys(ngrid*nlayer,zduadj)679 CALL zerophys(ngrid*nlayer,zdvadj)680 CALL zerophys(ngrid*nlayer,zdhadj)828 zduadj(:,:)=0 829 zdvadj(:,:)=0 830 zdhadj(:,:)=0 681 831 682 832 CALL convadj(ngrid,nlayer,nq,ptimestep, 683 $ pplay,pplev,zpopsk, 833 $ pplay,pplev,zpopsk,lmax_th, 684 834 $ pu,pv,zh,pq, 685 835 $ pdu,pdv,zdh,pdq, 686 836 $ zduadj,zdvadj,zdhadj, 687 837 $ zdqadj) 838 688 839 689 840 DO l=1,nlayer … … 718 869 $ co2ice,albedo,emis, 719 870 $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, 720 $ 871 $ fluxsurf_sw,zls) 721 872 722 873 DO l=1,nlayer … … 793 944 c ------------------ 794 945 946 #ifndef MESOSCALE 795 947 c -------------- 796 948 c photochemistry : … … 837 989 838 990 END IF ! of IF (photochem.or.thermochem) 991 #endif 839 992 840 993 c 7c. Aerosol particles … … 908 1061 endif ! of if (tracer) 909 1062 910 1063 #ifndef MESOSCALE 911 1064 c----------------------------------------------------------------------- 912 1065 c 8. THERMOSPHERE CALCULATION … … 933 1086 934 1087 endif ! of if (callthermos) 1088 #endif 935 1089 936 1090 c----------------------------------------------------------------------- … … 951 1105 952 1106 IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN 1107 #ifndef MESOSCALE 953 1108 if (caps.and.(obliquit.lt.27.)) then 954 1109 ! NB: Updated surface pressure, at grid point 'ngrid', is … … 957 1112 & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) 958 1113 endif 1114 #endif 959 1115 c ------------------------------------------------------------- 960 c Change of surface albedo (set to 0.4)in case of ground frost1116 c Change of surface albedo in case of ground frost 961 1117 c everywhere except on the north permanent cap and in regions 962 1118 c covered by dry ice. … … 965 1121 do ig=1,ngrid 966 1122 if ((co2ice(ig).eq.0).and. 967 & (qsurf(ig,igcm_h2o_ice).gt.0.005)) then 968 albedo(ig,1) = alb_surfice 969 albedo(ig,2) = alb_surfice 1123 & (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then 1124 albedo(ig,1) = albedo_h2o_ice 1125 albedo(ig,2) = albedo_h2o_ice 1126 c write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice) 1127 c write(*,*) "physiq.F frost :" 1128 c & ,lati(ig)*180./pi, long(ig)*180./pi 970 1129 endif 971 1130 enddo ! of do ig=1,ngrid … … 994 1153 995 1154 c----------------------------------------------------------------------- 996 c 1 1. Write output files1155 c 10. Write output files 997 1156 c ---------------------- 998 1157 … … 1039 1198 ENDDO 1040 1199 1041 c Compute surface stress : (NB: z0 is a common in planete.h) 1200 ! Potential Temperature 1201 1202 DO ig=1,ngridmx 1203 DO l=1,nlayermx 1204 zh(ig,l) = zt(ig,l)*(zplay(ig,l)/zplev(ig,1))**rcp 1205 ENDDO 1206 ENDDO 1207 1208 1209 c Compute surface stress : (NB: z0 is a common in surfdat.h) 1042 1210 c DO ig=1,ngrid 1043 c cd = (0.4/log(zzlay(ig,1)/z0 ))**21211 c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2 1044 1212 c zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2) 1045 1213 c ENDDO … … 1047 1215 c Sum of fluxes in solar spectral bands (for output only) 1048 1216 DO ig=1,ngrid 1049 1050 1217 fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2) 1218 fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2) 1051 1219 ENDDO 1052 1220 c ******* TEST ****************************************************** … … 1061 1229 ENDDO 1062 1230 ENDDO 1063 if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then 1231 if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then 1064 1232 write(*,*) 'PHYSIQ: stability WARNING :' 1065 1233 write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin), … … 1088 1256 1089 1257 IF (ngrid.NE.1) THEN 1090 print*,'Ls =',zls*180./pi, 1091 & ' tauref(700 Pa,lat=0) =',tauref(ngrid/2), 1092 & ' tau(Viking1) =',tau(ig_vl1,1) 1093 1094 1258 1259 #ifndef MESOSCALE 1095 1260 c ------------------------------------------------------------------- 1096 1261 c Writing NetCDF file "RESTARTFI" at the end of the run … … 1111 1276 . zgam,zthe) 1112 1277 ENDIF 1278 #endif 1113 1279 1114 1280 c ------------------------------------------------------------------- … … 1120 1286 if (water) then 1121 1287 1122 call zerophys(ngrid,mtot)1123 call zerophys(ngrid,icetot)1124 call zerophys(ngrid,rave)1125 call zerophys(ngrid,tauTES)1288 mtot(:)=0 1289 icetot(:)=0 1290 rave(:)=0 1291 tauTES(:)=0 1126 1292 do ig=1,ngrid 1127 1293 do l=1,nlayermx … … 1160 1326 c which can later be used to make the statistic files of the run: 1161 1327 c "stats") only possible in 3D runs ! 1162 1163 1328 1164 1329 IF (callstats) THEN … … 1268 1433 ENDIF 1269 1434 1435 1436 #ifdef MESOSCALE 1437 !!! 1438 !!! OUTPUT FIELDS 1439 !!! 1440 wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature 1441 wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice 1442 mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice 1443 !! JF 1444 TAU_lay(:)=tau(:,1)!!true opacity (not a reference like tauref) 1445 IF (igcm_dust_mass .ne. 0) THEN 1446 qsurfice_dust(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass) 1447 ENDIF 1448 IF (igcm_h2o_ice .ne. 0) THEN 1449 qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice) 1450 vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice) 1451 . * mugaz / mmol(igcm_h2o_ice) 1452 ENDIF 1453 !! Dust quantity integration along the vertical axe 1454 dustot(:)=0 1455 do ig=1,ngrid 1456 do l=1,nlayermx 1457 dustot(ig) = dustot(ig) + 1458 & zq(ig,l,igcm_dust_mass) 1459 & * (pplev(ig,l) - pplev(ig,l+1)) / g 1460 enddo 1461 enddo 1462 c AUTOMATICALLY GENERATED FROM REGISTRY 1463 #include "fill_save.inc" 1464 #else 1465 1270 1466 c ========================================================== 1271 1467 c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing … … 1275 1471 c WRITEDIAGFI can ALSO be called from any other subroutines 1276 1472 c for any variables !! 1277 call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, 1278 & emis) 1473 c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, 1474 c & emis) 1475 ! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay) 1476 ! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) 1279 1477 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, 1280 1478 & tsurf) … … 1282 1480 call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, 1283 1481 & co2ice) 1482 1483 c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7", 1484 c & "K",2,zt(1,7)) 1284 1485 c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, 1285 1486 c & fluxsurf_lw) … … 1291 1492 c & fluxtop_sw_tot) 1292 1493 call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) 1293 c call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)1294 1494 call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) 1295 1495 call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) 1296 1496 call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) 1297 ccall WRITEDIAGFI(ngrid,"rho","density","none",3,rho)1497 call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) 1298 1498 c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2) 1299 ccall WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)1499 ! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) 1300 1500 c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay) 1301 1501 c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, … … 1305 1505 c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', 1306 1506 c & 'w.m-2',3,zdtlw) 1307 1308 !!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL 1507 c CALL WRITEDIAGFI(ngridmx,'tauTESap', 1508 c & 'tau abs 825 cm-1', 1509 c & '',2,tauTES) 1510 1511 #ifdef MESOINI 1512 call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, 1513 & emis) 1514 call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) 1309 1515 call WRITEDIAGFI(ngrid,"tsoil","Soil temperature", 1310 1516 & "K",3,tsoil) 1311 1517 call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia", 1312 1518 & "K",3,inertiedat) 1313 !!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL 1519 #endif 1520 1314 1521 1315 1522 c ---------------------------------------------------------- … … 1324 1531 1325 1532 ! Compute co2 column 1326 c all zerophys(ngrid,co2col)1533 co2col(:)=0 1327 1534 do l=1,nlayermx 1328 1535 do ig=1,ngrid … … 1331 1538 enddo 1332 1539 enddo 1333 ccall WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,1334 c& co2col)1540 call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2, 1541 & co2col) 1335 1542 endif ! of if (tracer.and.(igcm_co2.ne.0)) 1336 1543 … … 1341 1548 if (water) then 1342 1549 1550 #ifdef MESOINI 1343 1551 !!!! waterice = q01, voir readmeteo.F90 1344 1552 call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice), … … 1353 1561 & 'kg.m-2',2, 1354 1562 & qsurf(1:ngridmx,igcm_h2o_ice)) 1563 #endif 1355 1564 1356 1565 CALL WRITEDIAGFI(ngridmx,'mtot', … … 1360 1569 & 'total mass of water ice', 1361 1570 & 'kg/m2',2,icetot) 1362 cc vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) 1363 cc & *mugaz/mmol(igcm_h2o_ice) 1364 cc call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr', 1365 cc & 'mol/mol',3,vmr) 1366 c CALL WRITEDIAGFI(ngridmx,'reffice', 1367 c & 'Mean reff', 1368 c & 'm',2,rave) 1369 cc call WRITEDIAGFI(ngridmx,'rice','Ice particle size', 1370 cc & 'm',3,rice) 1371 cc If activice is true, tauTES is computed in aeropacity.F; 1372 c if (.not.activice) then 1373 c CALL WRITEDIAGFI(ngridmx,'tauTESap', 1374 c & 'tau abs 825 cm-1', 1375 c & '',2,tauTES) 1376 c endif 1377 c call WRITEDIAGFI(ngridmx,'h2o_ice_s', 1378 c & 'surface h2o_ice', 1379 c & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) 1571 c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) 1572 c & *mugaz/mmol(igcm_h2o_ice) 1573 c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr', 1574 c & 'mol/mol',3,vmr) 1575 CALL WRITEDIAGFI(ngridmx,'reffice', 1576 & 'Mean reff', 1577 & 'm',2,rave) 1578 c call WRITEDIAGFI(ngridmx,'rice','Ice particle size', 1579 c & 'm',3,rice) 1580 c If activice is true, tauTES is computed in aeropacity.F; 1581 if (.not.activice) then 1582 CALL WRITEDIAGFI(ngridmx,'tauTESap', 1583 & 'tau abs 825 cm-1', 1584 & '',2,tauTES) 1585 endif 1586 call WRITEDIAGFI(ngridmx,'h2o_ice_s', 1587 & 'surface h2o_ice', 1588 & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) 1589 1590 if (caps) then 1591 do ig=1,ngridmx 1592 if (watercaptag(ig)) watercapflag(ig) = 1 1593 enddo 1594 CALL WRITEDIAGFI(ngridmx,'watercaptag', 1595 & 'Ice water caps', 1596 & '',2,watercapflag) 1597 endif 1598 CALL WRITEDIAGFI(ngridmx,'albedo', 1599 & 'albedo', 1600 & '',2,albedo(1:ngridmx,1)) 1380 1601 endif !(water) 1381 1602 … … 1407 1628 c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1)) 1408 1629 if (doubleq) then 1409 c ccall WRITEDIAGFI(ngridmx,'qsurf','qsurf',1410 c c& 'kg.m-2',2,qsurf(1,1))1411 c ccall WRITEDIAGFI(ngridmx,'Nsurf','N particles',1412 c c& 'N.m-2',2,qsurf(1,2))1413 c ccall WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',1414 c c& 'kg.m-2.s-1',2,zdqsdev(1,1))1415 c ccall WRITEDIAGFI(ngridmx,'dqssed','sedimentation',1416 c c& 'kg.m-2.s-1',2,zdqssed(1,1))1417 ccall WRITEDIAGFI(ngridmx,'reffdust','reffdust',1418 c& 'm',3,rdust*ref_r0)1630 c call WRITEDIAGFI(ngridmx,'qsurf','qsurf', 1631 c & 'kg.m-2',2,qsurf(1,1)) 1632 c call WRITEDIAGFI(ngridmx,'Nsurf','N particles', 1633 c & 'N.m-2',2,qsurf(1,2)) 1634 c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift', 1635 c & 'kg.m-2.s-1',2,zdqsdev(1,1)) 1636 c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation', 1637 c & 'kg.m-2.s-1',2,zdqssed(1,1)) 1638 call WRITEDIAGFI(ngridmx,'reffdust','reffdust', 1639 & 'm',3,rdust*ref_r0) 1419 1640 call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr', 1420 1641 & 'kg/kg',3,pq(1,1,igcm_dust_mass)) 1421 call WRITEDIAGFI(ngridmx,'dustN','Dust number', 1642 c call WRITEDIAGFI(ngridmx,'dustN','Dust number', 1643 c & 'part/kg',3,pq(1,1,igcm_dust_number)) 1644 #ifdef MESOINI 1645 call WRITEDIAGFI(ngridmx,'dustN','Dust number', 1422 1646 & 'part/kg',3,pq(1,1,igcm_dust_number)) 1647 #endif 1423 1648 else 1424 1649 do iq=1,dustbin … … 1437 1662 1438 1663 c ---------------------------------------------------------- 1664 c ---------------------------------------------------------- 1665 c PBL OUTPUS 1666 c ---------------------------------------------------------- 1667 c ---------------------------------------------------------- 1668 1669 1670 c ---------------------------------------------------------- 1671 c Outputs of surface layer 1672 c ---------------------------------------------------------- 1673 1674 1675 z_out=0. 1676 if (calltherm .and. (z_out .gt. 0.)) then 1677 call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th 1678 & ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar) 1679 1680 zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1)) 1681 call WRITEDIAGFI(ngridmx,'sqrt(zu2)', 1682 & 'horizontal velocity norm','m/s', 1683 & 2,zu2) 1684 1685 call WRITEDIAGFI(ngridmx,'Teta_out', 1686 & 'potential temperature at z_out','K', 1687 & 2,Teta_out) 1688 call WRITEDIAGFI(ngridmx,'u_out', 1689 & 'horizontal velocity norm at z_out','m/s', 1690 & 2,u_out) 1691 call WRITEDIAGFI(ngridmx,'u*', 1692 & 'friction velocity','m/s', 1693 & 2,ustar) 1694 call WRITEDIAGFI(ngridmx,'teta*', 1695 & 'friction potential temperature','K', 1696 & 2,tstar) 1697 else 1698 if((.not. calltherm).and.(z_out .gt. 0.)) then 1699 print*, 'WARNING : no interpolation in surface-layer :' 1700 print*, 'Outputing surface-layer quantities without thermals 1701 & does not make sense' 1702 endif 1703 endif 1704 1705 c ---------------------------------------------------------- 1706 c Outputs of thermals 1707 c ---------------------------------------------------------- 1708 if (calltherm) then 1709 1710 ! call WRITEDIAGFI(ngrid,'dtke', 1711 ! & 'tendance tke thermiques','m**2/s**2', 1712 ! & 3,dtke_th) 1713 ! call WRITEDIAGFI(ngrid,'d_u_ajs', 1714 ! & 'tendance u thermiques','m/s', 1715 ! & 3,pdu_th*ptimestep) 1716 ! call WRITEDIAGFI(ngrid,'d_v_ajs', 1717 ! & 'tendance v thermiques','m/s', 1718 ! & 3,pdv_th*ptimestep) 1719 ! if (tracer) then 1720 ! if (nq .eq. 2) then 1721 ! call WRITEDIAGFI(ngrid,'deltaq_th', 1722 ! & 'delta q thermiques','kg/kg', 1723 ! & 3,ptimestep*pdq_th(:,:,2)) 1724 ! endif 1725 ! endif 1726 1727 call WRITEDIAGFI(ngridmx,'lmax_th', 1728 & 'hauteur du thermique','K', 1729 & 2,lmax_th_out) 1730 call WRITEDIAGFI(ngridmx,'hfmax_th', 1731 & 'maximum TH heat flux','K.m/s', 1732 & 2,hfmax_th) 1733 call WRITEDIAGFI(ngridmx,'wmax_th', 1734 & 'maximum TH vertical velocity','m/s', 1735 & 2,wmax_th) 1736 1737 endif 1738 1739 c ---------------------------------------------------------- 1740 c ---------------------------------------------------------- 1741 c END OF PBL OUTPUS 1742 c ---------------------------------------------------------- 1743 c ---------------------------------------------------------- 1744 1745 1746 c ---------------------------------------------------------- 1439 1747 c Output in netcdf file "diagsoil.nc" for subterranean 1440 1748 c variables (output every "ecritphy", as for writediagfi) … … 1451 1759 c END OF WRITEDIAGFI 1452 1760 c ========================================================== 1761 #endif 1453 1762 1454 1763 ELSE ! if(ngrid.eq.1) … … 1470 1779 c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1') 1471 1780 1781 ! THERMALS STUFF 1D 1782 1783 z_out=0. 1784 if (calltherm .and. (z_out .gt. 0.)) then 1785 call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th 1786 & ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar) 1787 1788 zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1)) 1789 call WRITEDIAGFI(ngridmx,'sqrt(zu2)', 1790 & 'horizontal velocity norm','m/s', 1791 & 0,zu2) 1792 1793 call WRITEDIAGFI(ngridmx,'Teta_out', 1794 & 'potential temperature at z_out','K', 1795 & 0,Teta_out) 1796 call WRITEDIAGFI(ngridmx,'u_out', 1797 & 'horizontal velocity norm at z_out','m/s', 1798 & 0,u_out) 1799 call WRITEDIAGFI(ngridmx,'u*', 1800 & 'friction velocity','m/s', 1801 & 0,ustar) 1802 call WRITEDIAGFI(ngridmx,'teta*', 1803 & 'friction potential temperature','K', 1804 & 0,tstar) 1805 else 1806 if((.not. calltherm).and.(z_out .gt. 0.)) then 1807 print*, 'WARNING : no interpolation in surface-layer :' 1808 print*, 'Outputing surface-layer quantities without thermals 1809 & does not make sense' 1810 endif 1811 endif 1812 1813 if(calltherm) then 1814 1815 call WRITEDIAGFI(ngridmx,'lmax_th', 1816 & 'hauteur du thermique','point', 1817 & 0,lmax_th_out) 1818 call WRITEDIAGFI(ngridmx,'hfmax_th', 1819 & 'maximum TH heat flux','K.m/s', 1820 & 0,hfmax_th) 1821 call WRITEDIAGFI(ngridmx,'wmax_th', 1822 & 'maximum TH vertical velocity','m/s', 1823 & 0,wmax_th) 1824 1825 co2col(:)=0. 1826 if (tracer) then 1827 do l=1,nlayermx 1828 do ig=1,ngrid 1829 co2col(ig)=co2col(ig)+ 1830 & zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g 1831 enddo 1832 enddo 1833 1834 end if 1835 call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' & 1836 & ,'kg/m-2',0,co2col) 1837 endif 1838 call WRITEDIAGFI(ngrid,'w','vertical velocity' & 1839 & ,'m/s',1,pw) 1840 call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2) 1841 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0, 1842 & tsurf) 1843 call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu) 1844 call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv) 1845 1846 call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay) 1847 call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev) 1472 1848 ! or output in diagfi.nc (for testphys1d) 1473 1849 call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps) -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F
r609 r665 472 472 c-------------------------------------------------------------------- 473 473 474 justbackground=.false. 474 justbackground=.false. !! default value 475 write(*,*) "justbackground" 476 call getin("justbackground",justbackground) 477 write(*,*) "justbackground = ",justbackground 478 475 479 476 480 IF (justbackground .eq. .true.) THEN … … 775 779 DO l=1,nlayer 776 780 DO ig=1,ngrid 777 dsodust(ig,l) = dsodust(ig,l) + 778 & aerosol(ig,l,iaerdust(iaer)) * g / 779 & (pplev(ig,l) - pplev(ig,l+1)) 781 c dsodust(ig,l) = dsodust(ig,l) + 782 c & aerosol(ig,l,iaerdust(iaer)) * g / 783 c & (pplev(ig,l) - pplev(ig,l+1)) 784 dsodust(ig,l) = 785 & ( 0.75 * QREFvis3d(ig,l,iaer) / 786 & ( rho_dust * reffrad(ig,l,iaer) ) ) * 787 & pq(ig,l,igcm_dust_mass) 780 788 ENDDO 781 789 ENDDO -
trunk/MESOSCALE_DEV/MANUAL/SRC/faq.tex
r493 r665 195 195 %%% WPS PREP_MARS peuvent être liés entre e.g. pgf et mpi, ou ifort et mpifort 196 196 197 %%% BIG FILES declare -x WRFIO_NCD_LARGE_FILE_SUPPORT=1 198 197 199 \clearemptydoublepage -
trunk/UTIL/PYTHON/planetoplot.py
r647 r665 254 254 elif axtime == "sol": time = soltab 255 255 else: 256 if ope in ["cat"] and nnn > 0: count = time[-1] + 1 ## so that a cat is possible with simple subscripts 257 else: count = 0 256 258 if "Times" in nc.dimensions: time = count + np.arange(0,len(nc.dimensions["Times"]),1) 257 259 elif "Time" in nc.dimensions: time = count + np.arange(0,len(nc.dimensions["Time"]),1) 258 260 else: time = count + np.arange(0,1,1) 259 if ope in ["cat"]:260 if nnn > 0: count = time[-1] + 1 ## so that a cat is possible with simple subscripts261 else: count = 0262 else: count = 0263 261 if axtime in ["lt"]: 264 262 ftime = np.zeros(len(time))
Note: See TracChangeset
for help on using the changeset viewer.
