Changeset 56 for trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys
- Timestamp:
- Feb 13, 2011, 12:25:13 PM (14 years ago)
- Location:
- trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys
- Files:
-
- 1 added
- 1 edited
-
physiq.F (modified) (34 diffs)
-
physiq.F.newold (added)
Legend:
- Unmodified
- Added
- Removed
-
trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/physiq.F
r31 r56 30 30 c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff. 31 31 c 2. Compute radiative transfer tendencies 32 c (longwave and shortwave) for CO2 and dust.32 c (longwave and shortwave) for CO2 and aerosols. 33 33 c 3. Gravity wave and subgrid scale topography drag : 34 34 c 4. Vertical diffusion (turbulent mixing): … … 47 47 c - Dumping eof (if "calleofdump = .true.") 48 48 c - Output any needed variables in "diagfi" 49 c 10. Diagnostic: mass conservation of tracers49 c 11. Diagnostic: mass conservation of tracers 50 50 c 51 51 c author: … … 58 58 c Introduction of the thermosphere module: M. Angelats i Coll (2002) 59 59 c Water ice clouds: Franck Montmessin (update 06/2003) 60 c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009) 61 c Nb: See callradite.F for more information. 60 62 c 61 c62 c63 63 c arguments: 64 64 c ---------- … … 96 96 c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding 97 97 c pdt(ngrid,nlayermx) / variables due to physical processes. 98 c pdq(ngrid,nlayermx )/98 c pdq(ngrid,nlayermx,nqmx) / 99 99 c pdpsrf(ngrid) / 100 100 c tracerdyn call tracer in dynamical part of GCM ? … … 110 110 #include "comgeomfi.h" 111 111 #include "surfdat.h" 112 #include "comsoil.h" 112 113 #include "comdiurn.h" 113 114 #include "callkeys.h" … … 123 124 #include "chimiedata.h" 124 125 #include "watercap.h" 125 #include "fisice.h"126 126 #include "param.h" 127 127 #include "param_v3.h" … … 164 164 c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of 165 165 c aerosol optical properties : 166 REAL aerosol(ngridmx,nlayermx,naerkind) 166 REAL aerosol(ngridmx,nlayermx,naerkind) 167 167 168 168 INTEGER day_ini ! Initial date of the run (sol since Ls=0) … … 182 182 INTEGER ig_vl1 ! Grid Point near VL1 (for diagnostic) 183 183 184 c Variables used by the water ice microphysical scheme: 185 REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m) 186 REAL nuice(ngridmx,nlayermx) ! Estimated effective variance 187 ! of the size distribution 188 c Albedo of deposited surface ice 189 REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45 190 184 191 SAVE day_ini, icount 185 192 SAVE aerosol, tsurf,tsoil … … 195 202 c ----------------- 196 203 204 REAL CBRT 205 EXTERNAL CBRT 197 206 198 207 CHARACTER*80 fichier 199 INTEGER l,ig,ierr,igout,iq, tapphys 200 INTEGER iqmin ! Used if iceparty engaged 208 INTEGER l,ig,ierr,igout,iq,i, tapphys 201 209 202 210 REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2) … … 204 212 REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2) 205 213 REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2) 206 c for clear area (uncovered by clouds) :207 REAL clsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)208 REAL clsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)209 REAL cltop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)210 REAL cltop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)211 214 REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa 212 215 ! (used if active=F) … … 217 220 REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries 218 221 REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic) 219 220 222 221 223 c Tendancies due to various processes: … … 262 264 REAL ztime_fin 263 265 REAL zdh(ngridmx,nlayermx) 264 REAL pclc_min ! minimum of the cloud fraction over the whole domain265 266 INTEGER length 266 267 PARAMETER (length=100) … … 275 276 character*5 str5 276 277 real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx) 277 real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T) 278 REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei 279 ! (particules kg-1) 280 real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m) 278 281 real qtot1,qtot2 ! total aerosol mass 279 282 integer igmin, lmin 280 283 logical tdiag 281 284 282 285 real co2col(ngridmx) ! CO2 column 283 286 REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx) 284 287 REAL zstress(ngrid), cd 285 288 real hco2(nqmx),tmean, zlocal(nlayermx) 286 real rho(ngridmx,nlayermx) ! density 287 real vmr(ngridmx,nlayermx) ! volume mixing ratio 289 real rho(ngridmx,nlayermx) ! density 290 real vmr(ngridmx,nlayermx) ! volume mixing ratio 291 REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2) 292 REAL icetot(ngridmx) ! Total mass of water ice (kg/m2) 293 REAL rave(ngridmx) ! Mean water ice effective radius (m) 294 REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1 295 REAL tauTES(ngridmx) ! column optical depth at 825 cm-1 296 REAL Qabsice ! Water ice absorption coefficient 297 288 298 289 299 REAL time_phys 290 300 291 301 c======================================================================= 292 302 … … 334 344 s ptimestep,tsurf,tsoil,capcal,fluxgrd) 335 345 ELSE 336 PRINT*,'WARNING! Thermal conduction in the soil turned off' 346 PRINT*, 347 & 'PHYSIQ WARNING! Thermal conduction in the soil turned off' 337 348 DO ig=1,ngrid 338 349 capcal(ig)=1.e5 … … 378 389 enddo 379 390 endif 391 392 IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN 393 write(*,*)"physiq: water_param Surface ice alb:",alb_surfice 394 ENDIF 380 395 381 396 ENDIF ! (end of "if firstcall") … … 489 504 & CALL nlthermeq(ngrid,nlayer,pplev,pplay) 490 505 491 492 c Atmospheric dust opacity and aerosol distribution: 493 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 494 CALL dustopacity(ngrid,nlayer,nq,zday,pplay,pplev,zls,pq, 495 $ tauref,tau,aerosol) 506 c Note: Dustopacity.F has been transferred to callradite.F 496 507 497 508 c Call main radiative transfer scheme 498 509 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 499 c Transfer through dust and CO2, except NIR CO2 absorption 500 501 CALL callradite(icount,ngrid,nlayer,aerosol,albedo, 510 c Transfer through CO2 (except NIR CO2 absorption) 511 c and aerosols (dust and water ice) 512 513 c Radiative transfer 514 c ------------------ 515 CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, 502 516 $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, 503 $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw) 517 $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw, 518 & tauref,tau,aerosol,ccn,rdust,rice,nuice) 504 519 505 520 c CO2 near infrared absorption … … 538 553 ENDIF ! of if(mod(icount-1,iradia).eq.0) 539 554 540 541 542 555 c Transformation of the radiative tendencies: 543 556 c ------------------------------------------- … … 599 612 enddo 600 613 enddo 601 614 602 615 c Calling vdif (Martian version WITH CO2 condensation) 603 616 CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk, … … 700 713 $ co2ice,albedo,emis, 701 714 $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, 702 $ fluxsurf_sw )715 $ fluxsurf_sw,zls) 703 716 704 717 DO l=1,nlayer … … 738 751 c ---------------------------------------- 739 752 IF (water) THEN 740 call watercloud(ngrid,nlayer, ptimestep, 753 754 call watercloud(ngrid,nlayer,ptimestep, 741 755 & pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, 742 & pq,pdq,zdqcloud, qsurf,zdqscloud,zdtcloud,743 & nq,naerkind,tau, icount,zls)744 756 & pq,pdq,zdqcloud,zdqscloud,zdtcloud, 757 & nq,naerkind,tau, 758 & ccn,rdust,rice,nuice) 745 759 if (activice) then 746 760 c Temperature variation due to latent heat release … … 752 766 endif 753 767 754 IF (iceparty) then 755 iqmin=nq-1 756 ELSE 757 iqmin=nq 758 ENDIF 759 760 DO iq=iqmin,nq 768 ! increment water vapour and ice atmospheric tracers tendencies 769 IF (water) THEN 761 770 DO l=1,nlayer 762 DO ig=1,ngrid 763 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq) 764 ENDDO 771 DO ig=1,ngrid 772 pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+ 773 & zdqcloud(ig,l,igcm_h2o_vap) 774 pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+ 775 & zdqcloud(ig,l,igcm_h2o_ice) 776 ENDDO 765 777 ENDDO 766 DO ig=1,ngrid 767 dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq) 768 ENDDO 769 ENDDO 770 778 ENDIF ! of IF (water) THEN 779 ! Increment water ice surface tracer tendency 780 DO ig=1,ngrid 781 dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+ 782 & zdqscloud(ig,igcm_h2o_ice) 783 ENDDO 784 771 785 END IF ! of IF (water) 772 786 … … 779 793 IF (photochem .or. thermochem) then 780 794 call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, 781 . zzlay,zday,pq,pdq,zdqchim,zdqschim,zdqcloud,zdqscloud) 782 783 c Photochemistry includes condensation of H2O2 784 785 do iq=nqchem_min,nq 786 if (noms(iq).eq."h2o2") then 795 & zzlay,zday,pq,pdq,rice, 796 & zdqchim,zdqschim,zdqcloud,zdqscloud) 797 !NB: Photochemistry includes condensation of H2O2 798 799 ! increment values of tracers: 800 DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry 801 ! tracers is zero anyways 787 802 DO l=1,nlayer 788 803 DO ig=1,ngrid 789 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq) 790 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq) 804 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) 791 805 ENDDO 792 806 ENDDO 793 else 807 ENDDO ! of DO iq=1,nq 808 ! add condensation tendency for H2O2 809 if (igcm_h2o2.ne.0) then 794 810 DO l=1,nlayer 795 811 DO ig=1,ngrid 796 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq) 812 pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2) 813 & +zdqcloud(ig,l,igcm_h2o2) 797 814 ENDDO 798 815 ENDDO 799 endif 800 ENDDO 801 do iq=nqchem_min,nq 802 if (noms(iq).eq."h2o2") then 803 DO ig=1,ngrid 804 dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq) 805 dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq) 806 ENDDO 807 else 808 DO ig=1,ngrid 809 dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq) 810 ENDDO 811 endif 812 ENDDO 816 endif 817 818 ! increment surface values of tracers: 819 DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry 820 ! tracers is zero anyways 821 DO ig=1,ngrid 822 dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) 823 ENDDO 824 ENDDO ! of DO iq=1,nq 825 ! add condensation tendency for H2O2 826 if (igcm_h2o2.ne.0) then 827 DO ig=1,ngrid 828 dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2) 829 & +zdqscloud(ig,igcm_h2o2) 830 ENDDO 831 endif 813 832 814 833 END IF ! of IF (photochem.or.thermochem) … … 845 864 c ------------- 846 865 IF (sedimentation) THEN 847 call zerophys(ngrid*nlayer*nq, zdqsed) 848 call zerophys(ngrid*nq, zdqssed) 849 850 if(doubleq) then 851 call callsedim2q(ngrid,nlayer, ptimestep, 852 & pplev,zzlev, pt, 866 !call zerophys(ngrid*nlayer*nq, zdqsed) 867 zdqsed(1:ngrid,1:nlayer,1:nq)=0 868 !call zerophys(ngrid*nq, zdqssed) 869 zdqssed(1:ngrid,1:nq)=0 870 871 call callsedim(ngrid,nlayer, ptimestep, 872 & pplev,zzlev, pt, rdust, rice, 853 873 & pq, pdq, zdqsed, zdqssed,nq) 854 else855 call callsedim(ngrid,nlayer, ptimestep,856 & pplev,zzlev, pt,857 & pq, pdq, zdqsed, zdqssed,nq)858 end if859 860 874 861 875 DO iq=1, nq … … 899 913 & pt,pq,pu,pv,pdt,pdq, 900 914 $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff) 901 c do iq=nqchem_min,nq902 c write(*,*) 'thermo iq,pq',iq,pq(690,1,iq)903 c enddo904 915 905 916 DO l=1,nlayer … … 916 927 ENDDO 917 928 918 919 endif 929 endif ! of if (callthermos) 930 920 931 c----------------------------------------------------------------------- 921 932 c 9. Surface and sub-surface soil temperature … … 934 945 c corresponding to equilibrium temperature between phases of CO2 935 946 936 IF (tracer.AND.water.AND. ngridmx.NE.1) THEN947 IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN 937 948 if (caps.and.(obliquit.lt.27.)) then 938 949 ! NB: Updated surface pressure, at grid point 'ngrid', is … … 948 959 c ------------------------------------------------------------- 949 960 do ig=1,ngrid 950 951 c -------------- Version temporaire fit TES 2008 ------------ 952 if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then 953 albedo(ig,1)=0.45 954 albedo(ig,2)=0.45 955 endif 956 957 c if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then 958 c if (.not.watercaptag(ig)) then 959 c albedo(ig,1)=0.4 960 c albedo(ig,2)=0.4 961 c endif 962 c endif 963 c -------------- version Francois --------------- 964 c if (co2ice(ig).eq.0.and. 965 c & ((qsurf(ig,nqmx).gt.0.005).or.(watercaptag(ig)))) then 966 c albedo(ig,1)=max(0.4,albedodat(ig)) 967 c albedo(ig,2)=albedo(ig,1) 968 c endif 969 c --------------------------------------------- 961 if ((co2ice(ig).eq.0).and. 962 & (qsurf(ig,igcm_h2o_ice).gt.0.005)) then 963 albedo(ig,1) = alb_surfice 964 albedo(ig,2) = alb_surfice 965 endif 970 966 enddo ! of do ig=1,ngrid 971 ENDIF ! of IF (tracer.AND.water.AND. ngridmx.NE.1)967 ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1)) 972 968 973 969 c … … 1073 1069 ENDIF ! of IF (lwrite) 1074 1070 1075 1076 1071 IF (ngrid.NE.1) THEN 1077 1072 print*,'Ls =',zls*180./pi, 1078 1073 & ' tauref(700 Pa,lat=0) =',tauref(ngrid/2), 1079 1074 & ' tau(Viking1) =',tau(ig_vl1,1) 1075 1080 1076 1081 1077 c ------------------------------------------------------------------- … … 1098 1094 ENDIF 1099 1095 1096 c ------------------------------------------------------------------- 1097 c Calculation of diagnostic variables written in both stats and 1098 c diagfi files 1099 c ------------------------------------------------------------------- 1100 1101 if (tracer) then 1102 if (water) then 1103 1104 call zerophys(ngrid,mtot) 1105 call zerophys(ngrid,icetot) 1106 call zerophys(ngrid,rave) 1107 call zerophys(ngrid,tauTES) 1108 do ig=1,ngrid 1109 do l=1,nlayermx 1110 mtot(ig) = mtot(ig) + 1111 & zq(ig,l,igcm_h2o_vap) * 1112 & (pplev(ig,l) - pplev(ig,l+1)) / g 1113 icetot(ig) = icetot(ig) + 1114 & zq(ig,l,igcm_h2o_ice) * 1115 & (pplev(ig,l) - pplev(ig,l+1)) / g 1116 rave(ig) = rave(ig) + 1117 & zq(ig,l,igcm_h2o_ice) * 1118 & (pplev(ig,l) - pplev(ig,l+1)) / g * 1119 & rice(ig,l) * (1.+nuice_ref) 1120 c Computing abs optical depth at 825 cm-1 in each 1121 c layer to simulate NEW TES retrieval 1122 Qabsice = min( 1123 & max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2 1124 & ) 1125 opTES(ig,l)= 0.75 * Qabsice * 1126 & zq(ig,l,igcm_h2o_ice) * 1127 & (pplev(ig,l) - pplev(ig,l+1)) / g 1128 & / (rho_ice * rice(ig,l) * (1.+nuice_ref)) 1129 tauTES(ig)=tauTES(ig)+ opTES(ig,l) 1130 enddo 1131 rave(ig)=rave(ig)/max(icetot(ig),1.e-30) 1132 if (icetot(ig)*1e3.lt.0.01) rave(ig)=0. 1133 enddo 1134 1135 endif ! of if (water) 1136 endif ! of if (tracer) 1137 1100 1138 c ----------------------------------------------------------------- 1101 c Saving statistics :1139 c WSTATS: Saving statistics 1102 1140 c ----------------------------------------------------------------- 1103 1141 c ("stats" stores and accumulates 8 key variables in file "stats.nc" … … 1108 1146 IF (callstats) THEN 1109 1147 1110 1111 call wstats(ngrid,"ps","Surface pressure","K",2,ps) 1112 call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf) 1113 call wstats(ngrid,"co2ice","CO2 ice cover", 1114 . "kg.m-2",2,co2ice) 1115 c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, 1116 c . emis) 1117 call wstats(ngrid,"fluxsurf_lw", 1118 . "Thermal IR radiative flux to surface","W.m-2",2, 1119 . fluxsurf_lw) 1120 call wstats(ngrid,"fluxsurf_sw", 1121 . "Solar radiative flux to surface","W.m-2",2, 1122 . fluxsurf_sw_tot) 1123 call wstats(ngrid,"fluxtop_lw", 1124 . "Thermal IR radiative flux to space","W.m-2",2, 1125 . fluxtop_lw) 1126 call wstats(ngrid,"fluxtop_sw", 1127 . "Solar radiative flux to space","W.m-2",2, 1128 . fluxtop_sw_tot) 1129 call wstats(ngrid,"dod","Dust optical depth"," ",2,tau) 1130 1131 call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) 1132 call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) 1133 call wstats(ngrid,"v","Meridional (North-South) wind", 1134 . "m.s-1",3,zv) 1135 call wstats(ngrid,"w","Vertical (down-up) wind", 1136 . "m.s-1",3,pw) 1137 call wstats(ngrid,"rho","Atmospheric density","none",3,rho) 1138 call wstats(ngrid,"q2", 1139 . "Boundary layer eddy kinetic energy","m2.s-2",3,q2) 1148 call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) 1149 call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf) 1150 call wstats(ngrid,"co2ice","CO2 ice cover", 1151 & "kg.m-2",2,co2ice) 1152 call wstats(ngrid,"fluxsurf_lw", 1153 & "Thermal IR radiative flux to surface","W.m-2",2, 1154 & fluxsurf_lw) 1155 call wstats(ngrid,"fluxsurf_sw", 1156 & "Solar radiative flux to surface","W.m-2",2, 1157 & fluxsurf_sw_tot) 1158 call wstats(ngrid,"fluxtop_lw", 1159 & "Thermal IR radiative flux to space","W.m-2",2, 1160 & fluxtop_lw) 1161 call wstats(ngrid,"fluxtop_sw", 1162 & "Solar radiative flux to space","W.m-2",2, 1163 & fluxtop_sw_tot) 1164 call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) 1165 call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) 1166 call wstats(ngrid,"v","Meridional (North-South) wind", 1167 & "m.s-1",3,zv) 1168 call wstats(ngrid,"w","Vertical (down-up) wind", 1169 & "m.s-1",3,pw) 1170 call wstats(ngrid,"rho","Atmospheric density","none",3,rho) 1171 call wstats(ngrid,"pressure","Pressure","Pa",3,pplay) 1172 c call wstats(ngrid,"q2", 1173 c & "Boundary layer eddy kinetic energy", 1174 c & "m2.s-2",3,q2) 1175 c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, 1176 c & emis) 1177 c call wstats(ngrid,"ssurf","Surface stress","N.m-2", 1178 c & 2,zstress) 1179 c call wstats(ngrid,"sw_htrt","sw heat.rate", 1180 c & "W.m-2",3,zdtsw) 1181 c call wstats(ngrid,"lw_htrt","lw heat.rate", 1182 c & "W.m-2",3,zdtlw) 1140 1183 1141 1184 if (tracer) then 1142 if (water) then 1143 vmr=zq(1:ngridmx,1:nlayermx,nqmx)*mugaz/mmol(nqmx) 1185 if (water) then 1186 vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap) 1187 & *mugaz/mmol(igcm_h2o_vap) 1144 1188 call wstats(ngrid,"vmr_h2ovapor", 1145 . "H2O vapor volume mixing ratio","mol/mol", 1146 . 3,vmr) 1147 if (iceparty) then 1148 vmr=zq(1:ngridmx,1:nlayermx,nqmx-1)*mugaz/mmol(nqmx-1) 1149 call wstats(ngrid,"vmr_h2oice", 1150 . "H2O ice volume mixing ratio","mol/mol", 1151 . 3,vmr) 1189 & "H2O vapor volume mixing ratio","mol/mol", 1190 & 3,vmr) 1191 vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) 1192 & *mugaz/mmol(igcm_h2o_ice) 1193 call wstats(ngrid,"vmr_h2oice", 1194 & "H2O ice volume mixing ratio","mol/mol", 1195 & 3,vmr) 1196 call wstats(ngrid,"h2o_ice_s", 1197 & "surface h2o_ice","kg/m2", 1198 & 2,qsurf(1,igcm_h2o_ice)) 1199 1200 call wstats(ngrid,"mtot", 1201 & "total mass of water vapor","kg/m2", 1202 & 2,mtot) 1203 call wstats(ngrid,"icetot", 1204 & "total mass of water ice","kg/m2", 1205 & 2,icetot) 1206 call wstats(ngrid,"reffice", 1207 & "Mean reff","m", 1208 & 2,rave) 1209 c call wstats(ngrid,"rice", 1210 c & "Ice particle size","m", 1211 c & 3,rice) 1212 c If activice is true, tauTES is computed in aeropacity.F; 1213 if (.not.activice) then 1214 call wstats(ngrid,"tauTESap", 1215 & "tau abs 825 cm-1","", 1216 & 2,tauTES) 1152 1217 endif 1153 endif 1218 1219 endif ! of if (water) 1154 1220 1155 1221 if (thermochem.or.photochem) then 1156 1222 do iq=1,nq 1157 1223 if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or. 1158 . (noms(iq).eq."co").or.(noms(iq).eq."n2").or.1159 . (noms(iq).eq."h2").or.1160 . (noms(iq).eq."o3")) then1224 . (noms(iq).eq."co").or.(noms(iq).eq."n2").or. 1225 . (noms(iq).eq."h2").or. 1226 . (noms(iq).eq."o3")) then 1161 1227 do l=1,nlayer 1162 1228 do ig=1,ngrid … … 1168 1234 endif 1169 1235 enddo 1170 endif 1171 endif !tracer 1172 1173 IF(lastcall) THEN 1174 write (*,*) "Writing stats..." 1175 call mkstats(ierr) 1176 ENDIF 1177 ENDIF !if callstats 1236 endif ! of if (thermochem.or.photochem) 1237 1238 endif ! of if (tracer) 1239 1240 IF(lastcall) THEN 1241 write (*,*) "Writing stats..." 1242 call mkstats(ierr) 1243 ENDIF 1244 1245 ENDIF !if callstats 1178 1246 1179 1247 c (Store EOF for Mars Climate database software) … … 1182 1250 ENDIF 1183 1251 1184 1185 c ----------------------------------------------------------------------1186 c output in netcdf file "DIAGFI", containing any variable for diagnostic1187 c (output with period"ecritphy", set in "run.def")1188 c ----------------------------------------------------------------------1252 c ========================================================== 1253 c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing 1254 c any variable for diagnostic (output with period 1255 c "ecritphy", set in "run.def") 1256 c ========================================================== 1189 1257 c WRITEDIAGFI can ALSO be called from any other subroutines 1190 1258 c for any variables !! 1191 call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,1192 .emis)1193 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,1194 .tsurf)1195 call WRITEDIAGFI(ngrid,"ps","surface pressure","K",2,ps)1259 call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, 1260 & emis) 1261 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, 1262 & tsurf) 1263 call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) 1196 1264 call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, 1197 .co2ice)1198 c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,1199 c .fluxsurf_lw)1200 c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,1201 c .fluxsurf_sw_tot)1202 c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,1203 c .fluxtop_lw)1204 c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,1205 c .fluxtop_sw_tot)1206 call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)1265 & co2ice) 1266 c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, 1267 c & fluxsurf_lw) 1268 c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2, 1269 c & fluxsurf_sw_tot) 1270 c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2, 1271 c & fluxtop_lw) 1272 c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2, 1273 c & fluxtop_sw_tot) 1274 call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) 1207 1275 c call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau) 1208 1276 call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) 1209 1277 call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) 1210 1278 call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) 1211 c call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)1279 c call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) 1212 1280 c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2) 1213 c call WRITEDIAGFI(ngridm,'Teta','T potentielle','K',3,zh) 1214 c call WRITEDIAGFI(ngridm,'Pression','Pression','Pa',3,pplay) 1215 call WRITEDIAGFI(ngrid,"tsoil","soil temperature","K",3,tsoil) 1216 1217 1218 c OUTPUT of tracer mass mixing ratio and surface value : 1219 if (tracer) then 1220 c (for photochemistry, outputs are done in calchim) 1221 do iq=1,nqmx 1222 write(str2(1:2),'(i2.2)') iq 1223 call WRITEDIAGFI(ngridmx,'q'//str2,noms(iq), 1281 c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) 1282 c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay) 1283 c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, 1284 c & zstress) 1285 c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate', 1286 c & 'w.m-2',3,zdtsw) 1287 c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', 1288 c & 'w.m-2',3,zdtlw) 1289 1290 !!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL 1291 call WRITEDIAGFI(ngrid,"tsoil","Soil temperature", 1292 & "K",3,tsoil) 1293 call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia", 1294 & "K",3,inertiedat) 1295 !!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL 1296 1297 c ---------------------------------------------------------- 1298 c Outputs of the CO2 cycle 1299 c ---------------------------------------------------------- 1300 1301 if (tracer.and.(igcm_co2.ne.0)) then 1302 ! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer", 1303 ! & "kg/kg",2,zq(1,1,igcm_co2)) 1304 ! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio", 1305 ! & "kg/kg",3,zq(1,1,igcm_co2)) 1306 1307 ! Compute co2 column 1308 call zerophys(ngrid,co2col) 1309 do l=1,nlayermx 1310 do ig=1,ngrid 1311 co2col(ig)=co2col(ig)+ 1312 & zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g 1313 enddo 1314 enddo 1315 c call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2, 1316 c & co2col) 1317 endif ! of if (tracer.and.(igcm_co2.ne.0)) 1318 1319 c ---------------------------------------------------------- 1320 c Outputs of the water cycle 1321 c ---------------------------------------------------------- 1322 if (tracer) then 1323 if (water) then 1324 1325 !!!! waterice = q01, voir readmeteo.F90 1326 call WRITEDIAGFI(ngridmx,'q01',noms(iq), 1327 & 'kg/kg',3, 1328 & zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)) 1329 !!!! watervapor = q02, voir readmeteo.F90 1330 call WRITEDIAGFI(ngridmx,'q02',noms(iq), 1331 & 'kg/kg',3, 1332 & zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)) 1333 1334 c call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq), 1335 c & 'kg.m-2',2,qsurf(1,iq)) 1336 1337 1338 c CALL WRITEDIAGFI(ngridmx,'mtot', 1339 c & 'total mass of water vapor', 1340 c & 'kg/m2',2,mtot) 1341 c CALL WRITEDIAGFI(ngridmx,'icetot', 1342 c & 'total mass of water ice', 1343 c & 'kg/m2',2,icetot) 1344 cc vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) 1345 cc & *mugaz/mmol(igcm_h2o_ice) 1346 cc call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr', 1347 cc & 'mol/mol',3,vmr) 1348 c CALL WRITEDIAGFI(ngridmx,'reffice', 1349 c & 'Mean reff', 1350 c & 'm',2,rave) 1351 cc call WRITEDIAGFI(ngridmx,'rice','Ice particle size', 1352 cc & 'm',3,rice) 1353 cc If activice is true, tauTES is computed in aeropacity.F; 1354 c if (.not.activice) then 1355 c CALL WRITEDIAGFI(ngridmx,'tauTESap', 1356 c & 'tau abs 825 cm-1', 1357 c & '',2,tauTES) 1358 c endif 1359 c call WRITEDIAGFI(ngridmx,'h2o_ice_s', 1360 c & 'surface h2o_ice', 1361 c & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) 1362 endif !(water) 1363 1364 1365 if (water.and..not.photochem) then 1366 iq=nq 1367 c write(str2(1:2),'(i2.2)') iq 1368 c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud', 1369 c & 'kg.m-2',2,zdqscloud(1,iq)) 1370 c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim', 1371 c & 'kg/kg',3,zdqchim(1,1,iq)) 1372 c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif', 1373 c & 'kg/kg',3,zdqdif(1,1,iq)) 1374 c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj', 1375 c & 'kg/kg',3,zdqadj(1,1,iq)) 1376 c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c', 1377 c & 'kg/kg',3,zdqc(1,1,iq)) 1378 endif !(water.and..not.photochem) 1379 endif 1380 1381 c ---------------------------------------------------------- 1382 c Outputs of the dust cycle 1383 c ---------------------------------------------------------- 1384 1385 c call WRITEDIAGFI(ngridmx,'tauref', 1386 c & 'Dust ref opt depth','NU',2,tauref) 1387 1388 if (tracer.and.(dustbin.ne.0)) then 1389 c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1)) 1390 if (doubleq) then 1391 cc call WRITEDIAGFI(ngridmx,'qsurf','qsurf', 1392 cc & 'kg.m-2',2,qsurf(1,1)) 1393 cc call WRITEDIAGFI(ngridmx,'Nsurf','N particles', 1394 cc & 'N.m-2',2,qsurf(1,2)) 1395 cc call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift', 1396 cc & 'kg.m-2.s-1',2,zdqsdev(1,1)) 1397 cc call WRITEDIAGFI(ngridmx,'dqssed','sedimentation', 1398 cc & 'kg.m-2.s-1',2,zdqssed(1,1)) 1399 c call WRITEDIAGFI(ngridmx,'reffdust','reffdust', 1400 c & 'm',3,rdust*ref_r0) 1401 call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr', 1402 & 'kg/kg',3,pq(1,1,igcm_dust_mass)) 1403 call WRITEDIAGFI(ngridmx,'dustN','Dust number', 1404 & 'part/kg',3,pq(1,1,igcm_dust_number)) 1405 else 1406 do iq=1,dustbin 1407 write(str2(1:2),'(i2.2)') iq 1408 call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio', 1224 1409 & 'kg/kg',3,zq(1,1,iq)) 1225 call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq), 1226 & 'kg.m-2',2,qsurf(1,iq)) 1227 end do 1228 end if 1229 c *************************************************** 1230 1231 c Outputs for H2O 1232 if (tracer) then 1233 if (activice) then 1234 c call WRITEDIAGFI(ngridmx,'tauice','tau','SI',2,tau(1,2)) 1235 c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate', 1236 c & 'w.m-2',3,zdtsw) 1237 c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', 1238 c & 'w.m-2',3,zdtlw) 1239 endif !(activice) 1240 1241 if (water.and..not.photochem) then 1242 iq=nq 1243 c write(str2(1:2),'(i2.2)') iq 1244 c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud', 1245 c & 'kg.m-2',2,zdqscloud(1,iq)) 1246 c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim', 1247 c & 'kg/kg',3,zdqchim(1,1,iq)) 1248 c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif', 1249 c & 'kg/kg',3,zdqdif(1,1,iq)) 1250 c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj', 1251 c & 'kg/kg',3,zdqadj(1,1,iq)) 1252 c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c', 1253 c & 'kg/kg',3,zdqc(1,1,iq)) 1254 endif !(water.and..not.photochem) 1255 1256 c if (iceparty) then 1257 c iq=nq-1 1258 c write(str2(1:2),'(i2.2)') iq 1259 c call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio', 1260 c & 'kg/kg',3,zq(1,1,iq)) 1261 c endif !(iceparty) 1262 endif 1263 1264 c Outputs for dust tracers 1265 1266 if (tracer.and.(dustbin.ne.0)) then 1267 call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1)) 1268 if (doubleq) then 1269 call WRITEDIAGFI(ngridmx,'qsurf','qsurf', 1270 & 'kg.m-2',2,qsurf(1,1)) 1271 call WRITEDIAGFI(ngridmx,'Nsurf','N particles', 1272 & 'N.m-2',2,qsurf(1,2)) 1273 call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift', 1274 & 'kg.m-2.s-1',2,zdqsdev(1,1)) 1275 call WRITEDIAGFI(ngridmx,'dqssed','sedimentation', 1276 & 'kg.m-2.s-1',2,zdqssed(1,1)) 1277 do l=1,nlayer 1278 do ig=1, ngrid 1279 reff(ig,l)= ref_r0 * 1280 & (r3n_q*pq(ig,l,1)/max(pq(ig,l,2),0.01))**(1./3.) 1281 reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6) 1282 end do 1283 end do 1284 call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff) 1285 else 1286 do iq=1,dustbin 1287 write(str2(1:2),'(i2.2)') iq 1288 call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio', 1289 & 'kg/kg',3,zq(1,1,iq)) 1290 call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf', 1291 & 'kg.m-2',2,qsurf(1,iq)) 1292 end do 1293 endif ! (doubleq) 1410 call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf', 1411 & 'kg.m-2',2,qsurf(1,iq)) 1412 end do 1413 endif ! (doubleq) 1414 c if (submicron) then 1415 c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr', 1416 c & 'kg/kg',3,pq(1,1,igcm_dust_submicron)) 1417 c endif ! (submicron) 1294 1418 end if ! (tracer.and.(dustbin.ne.0)) 1295 1419 1420 c ---------------------------------------------------------- 1421 c Output in netcdf file "diagsoil.nc" for subterranean 1422 c variables (output every "ecritphy", as for writediagfi) 1423 c ---------------------------------------------------------- 1424 1425 ! Write soil temperature 1426 ! call writediagsoil(ngrid,"soiltemp","Soil temperature","K", 1427 ! & 3,tsoil) 1428 ! Write surface temperature 1429 ! call writediagsoil(ngrid,"tsurf","Surface temperature","K", 1430 ! & 2,tsurf) 1431 1432 c ========================================================== 1433 c END OF WRITEDIAGFI 1434 c ========================================================== 1296 1435 1297 1436 ELSE ! if(ngrid.eq.1) 1298 1437 1299 1438 print*,'Ls =',zls*180./pi, 1300 & ' tauref(700 Pa) =',tauref ,' local tau =',tau(1,1)1439 & ' tauref(700 Pa) =',tauref 1301 1440 c ---------------------------------------------------------------------- 1302 1441 c Output in grads file "g1d" (ONLY when using testphys1d) … … 1305 1444 c 1306 1445 c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2') 1307 c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')1308 CALL writeg1d(ngrid,1,ps,'ps','Pa')1446 c CALL writeg1d(ngrid,1,tsurf,'tsurf','K') 1447 c CALL writeg1d(ngrid,1,ps,'ps','Pa') 1309 1448 1310 CALL writeg1d(ngrid,nlayer,zt,'T','K')1449 c CALL writeg1d(ngrid,nlayer,zt,'T','K') 1311 1450 c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1') 1312 1451 c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1') 1313 1452 c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1') 1314 1453 1454 ! or output in diagfi.nc (for testphys1d) 1455 call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps) 1456 call WRITEDIAGFI(ngridmx,'temp','Temperature', 1457 & 'K',1,zt) 1458 1315 1459 if(tracer) then 1316 1460 c CALL writeg1d(ngrid,1,tau,'tau','SI') 1317 1461 do iq=1,nq 1318 CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 1462 c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 1463 call WRITEDIAGFI(ngridmx,trim(noms(iq)), 1464 & trim(noms(iq)),'kg/kg',1,zq(1,1,iq)) 1319 1465 end do 1320 1466 end if … … 1333 1479 & rnew(1,nlayer)*tmean/g 1334 1480 1335 c if(tracer) then1336 c do l=2,nlayer1337 c do iq=1,nq1338 c hco2(iq)=zq(1,l,iq)/zq(1,l-1,iq)1339 c hco2(iq)=-(zlocal(l)-zlocal(l-1))/log(hco2(iq))/1000.1340 c enddo1341 c write(225,*) l,pt(1,l),(hco2(iq),iq=1,6)1342 c write(226,*) l,(hco2(iq),iq=7,13)1343 c enddo1344 c endif1345 1346 1481 END IF ! if(ngrid.ne.1) 1347 1348 1482 1349 1483 icount=icount+1
Note: See TracChangeset
for help on using the changeset viewer.
