Changeset 2931 for trunk/LMDZ.MARS
- Timestamp:
- Apr 3, 2023, 4:49:33 PM (20 months ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 3 edited
-
README (modified) (1 diff)
-
libf/phymars/physiq_mod.F (modified) (16 diffs)
-
libf/phymars/topmons_mod.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r2930 r2931 3979 3979 == 02/04/2023 == EM 3980 3980 Adapt topmons_setup to handle both cartesian and icosahedral grids 3981 3982 ==03/04/2023 == RV 3983 Remove the test ngrid==1 for writediagfi. The routine can handle well the 1D case. -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r2926 r2931 2512 2512 2513 2513 #ifndef MESOSCALE 2514 IF (ngrid .eq. 1) THEN2515 dimout=02516 ELSE2517 dimout=22518 ENDIF2519 2514 DO n=1,n_out 2520 2515 write(zstring, '(F8.6)') z_out(n) 2521 2516 call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring), 2522 & 'potential temperature at z_out','K', dimout,T_out(:,n))2517 & 'potential temperature at z_out','K',2,T_out(:,n)) 2523 2518 call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring), 2524 & 'horizontal velocity norm at z_out','m/s', dimout,u_out(:,n))2519 & 'horizontal velocity norm at z_out','m/s',2,u_out(:,n)) 2525 2520 ENDDO 2526 2521 call WRITEDIAGFI(ngrid,'u_star', 2527 & 'friction velocity','m/s', dimout,ustar)2522 & 'friction velocity','m/s',2,ustar) 2528 2523 call WRITEDIAGFI(ngrid,'teta_star', 2529 & 'friction potential temperature','K',dimout,tstar) 2530 ! call WRITEDIAGFI(ngrid,'L', 2531 ! & 'Monin Obukhov length','m',dimout,L_mo) 2524 & 'friction potential temperature','K',2,tstar) 2532 2525 call WRITEDIAGFI(ngrid,'vvv', 2533 & 'Vertical velocity variance at zout','m', dimout,vvv)2526 & 'Vertical velocity variance at zout','m',2,vvv) 2534 2527 call WRITEDIAGFI(ngrid,'vhf', 2535 & 'Vertical heat flux at zout','m', dimout,vhf)2528 & 'Vertical heat flux at zout','m',2,vhf) 2536 2529 #else 2537 2530 T_out1(:)=T_out(:,1) … … 2610 2603 #endif 2611 2604 2612 IF (ngrid.NE.1) then2605 c IF (ngrid.NE.1) then 2613 2606 2614 2607 c ------------------------------------------------------------------- … … 2640 2633 if(doubleq) then 2641 2634 do ig=1,ngrid 2635 IF (sedimentation) THEN 2642 2636 dqdustsurf(ig) = 2643 2637 & zdqssed(ig,igcm_dust_mass)*tauscaling(ig) 2644 2638 dndustsurf(ig) = 2645 2639 & zdqssed(ig,igcm_dust_number)*tauscaling(ig) 2640 ENDIF 2646 2641 ndust(ig,:) = 2647 2642 & zq(ig,:,igcm_dust_number)*tauscaling(ig) … … 2651 2646 if (scavenging) then 2652 2647 do ig=1,ngrid 2648 IF (sedimentation) THEN 2653 2649 dqdustsurf(ig) = dqdustsurf(ig) + 2654 2650 & zdqssed(ig,igcm_ccn_mass)*tauscaling(ig) 2655 2651 dndustsurf(ig) = dndustsurf(ig) + 2656 2652 & zdqssed(ig,igcm_ccn_number)*tauscaling(ig) 2653 ENDIF 2657 2654 nccn(ig,:) = 2658 2655 & zq(ig,:,igcm_ccn_number)*tauscaling(ig) … … 2853 2850 call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, 2854 2851 & emis(:,iflat)) 2855 c call wstats(ngrid,"ssurf","Surface stress","N.m-2",2856 c & 2,zstress)2857 c call wstats(ngrid,"sw_htrt","sw heat.rate",2858 c & "W.m-2",3,zdtsw)2859 c call wstats(ngrid,"lw_htrt","lw heat.rate",2860 c & "W.m-2",3,zdtlw)2861 2852 call wstats(ngrid,"fluxsurf_dir_dn_sw", 2862 2853 & "Direct incoming SW flux at surface", … … 2950 2941 2951 2942 if (doubleq) then 2952 c call wstats(ngrid,'qsurf','qsurf',2953 c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))2954 c call wstats(ngrid,'Nsurf','N particles',2955 c & 'N.m-2',2,qsurf(1,igcm_dust_number))2956 c call wstats(ngrid,'dqsdev','ddevil lift',2957 c & 'kg.m-2.s-1',2,zdqsdev(1,1))2958 c call wstats(ngrid,'dqssed','sedimentation',2959 c & 'kg.m-2.s-1',2,zdqssed(1,1))2960 c call wstats(ngrid,'dqsdif','diffusion',2961 c & 'kg.m-2.s-1',2,zdqsdif(1,1))2962 2943 call wstats(ngrid,'dqsdust', 2963 2944 & 'deposited surface dust mass', … … 3159 3140 & "Surface emissivity","w.m-1",2, emis(:,islope)) 3160 3141 ENDDO 3161 c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)3162 c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)3163 3142 call WRITEDIAGFI(ngrid,"zzlev","Interlayer altitude", 3164 3143 & "m",3,zzlev(:,1:nlayer)) … … 3223 3202 call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) 3224 3203 call WRITEDIAGFI(ngrid,"rho","density","kg.m-3",3,rho) 3225 c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)3226 c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)3227 3204 call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay) 3228 c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, 3229 c & zstress) 3205 call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) 3230 3206 call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate', 3231 3207 & 'K/s',3,zdtsw) … … 3501 3477 end if !(water) 3502 3478 3503 if (water.and..not.photochem) then3504 iq = nq3505 c write(str2(1:2),'(i2.2)') iq3506 c call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',3507 c & 'kg.m-2',2,zdqscloud(1,iq))3508 c call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',3509 c & 'kg/kg',3,zdqchim(1,1,iq))3510 c call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',3511 c & 'kg/kg',3,zdqdif(1,1,iq))3512 c call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',3513 c & 'kg/kg',3,zdqadj(1,1,iq))3514 c call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',3515 c & 'kg/kg',3,zdqc(1,1,iq))3516 end if !(water.and..not.photochem)3517 3518 3479 c ---------------------------------------------------------- 3519 3480 c Outputs of the dust cycle … … 3540 3501 #ifndef MESOINI 3541 3502 if (doubleq) then 3542 c call WRITEDIAGFI(ngrid,'qsurf','qsurf',3543 c & 'kg.m-2',2,qsurf(1,igcm_dust_mass))3544 c call WRITEDIAGFI(ngrid,'Nsurf','N particles',3545 c & 'N.m-2',2,qsurf(1,igcm_dust_number))3546 c call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',3547 c & 'kg.m-2.s-1',2,zdqsdev(1,1))3548 c call WRITEDIAGFI(ngrid,'dqssed','sedimentation',3549 c & 'kg.m-2.s-1',2,zdqssed(1,1))3550 c call WRITEDIAGFI(ngrid,'dqsdif','diffusion',3551 c & 'kg.m-2.s-1',2,zdqsdif(1,1))3552 c call WRITEDIAGFI(ngrid,'sedice','sedimented ice',3553 c & 'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))3554 c call WRITEDIAGFI(ngrid,'subice','sublimated ice',3555 c & 'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))3556 3503 call WRITEDIAGFI(ngrid,'dqsdust', 3557 3504 & 'deposited surface dust mass', … … 3648 3595 & ,'sedimentation q','kg.m-2.s-1',3, 3649 3596 & zdqsed(:,:,igcm_stormdust_mass)) 3597 call WRITEDIAGFI(ngrid,'zdqsed_dustN' 3598 & ,'sedimentation N','Nbr.m-2.s-1',3, 3599 & zdqsed(:,:,igcm_dust_number)) 3650 3600 call WRITEDIAGFI(ngrid,'rdust','rdust', 3651 3601 & 'm',3,rdust) … … 3714 3664 ENDDO 3715 3665 endif ! (scavenging) 3716 3717 c if (submicron) then3718 c call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',3719 c & 'kg/kg',3,pq(1,1,igcm_dust_submicron))3720 c endif ! (submicron)3721 3666 3722 3667 #else … … 3789 3734 $ 0,[PhiEscD]) 3790 3735 3791 ! call wstats(ngrid,"PhiH","H escape flux","s-1",3792 ! $ 0,[PhiEscH])3793 ! call wstats(ngrid,"PhiH2","H2 escape flux","s-1",3794 ! $ 0,[PhiEscH2])3795 ! call wstats(ngrid,"PhiD","D escape flux","s-1",3796 ! $ 0,[PhiEscD])3797 3798 ! call wstats(ngrid,"q15um","15 um cooling","K/s",3799 ! $ 3,zdtnlte)3800 ! call wstats(ngrid,"quv","UV heating","K/s",3801 ! $ 3,zdteuv)3802 ! call wstats(ngrid,"cond","Thermal conduction","K/s",3803 ! $ 3,zdtconduc)3804 ! call wstats(ngrid,"qnir","NIR heating","K/s",3805 ! $ 3,zdtnirco2)3806 3807 3736 endif !(callthermos) 3808 3737 … … 3822 3751 c ---------------------------------------------------------- 3823 3752 if (calltherm) then 3824 ! call WRITEDIAGFI(ngrid,'dtke', 3825 ! & 'tendance tke thermiques','m**2/s**2', 3826 ! & 3,dtke_th) 3827 ! call WRITEDIAGFI(ngrid,'d_u_ajs', 3828 ! & 'tendance u thermiques','m/s', 3829 ! & 3,pdu_th*ptimestep) 3830 ! call WRITEDIAGFI(ngrid,'d_v_ajs', 3831 ! & 'tendance v thermiques','m/s', 3832 ! & 3,pdv_th*ptimestep) 3833 ! if (nq .eq. 2) then 3834 ! call WRITEDIAGFI(ngrid,'deltaq_th', 3835 ! & 'delta q thermiques','kg/kg', 3836 ! & 3,ptimestep*pdq_th(:,:,2)) 3837 ! end if 3838 3753 call WRITEDIAGFI(ngrid,'lmax_th', 3754 & 'hauteur du thermique','point', 3755 & 2,lmax_th_out) 3839 3756 call WRITEDIAGFI(ngrid,'zmax_th', 3840 3757 & 'hauteur du thermique','m', … … 3869 3786 & 3,tsoil(:,:,islope)) 3870 3787 ENDDO 3871 ! Write surface temperature 3872 ! call writediagsoil(ngrid,"tsurf","Surface temperature","K", 3873 ! & 2,tsurf) 3874 3875 c ========================================================== 3876 c END OF WRITEDIAGFI 3877 c ========================================================== 3878 #endif 3879 ! of ifdef MESOSCALE 3880 3881 ELSE ! if(ngrid.eq.1) 3882 3883 #ifndef MESOSCALE 3884 write(*, 3885 & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)') 3886 & zls*180./pi,odpref,tau_pref_scenario 3887 c ---------------------------------------------------------------------- 3888 c Output in grads file "g1d" (ONLY when using testphys1d) 3889 c (output at every X physical timestep) 3890 c ---------------------------------------------------------------------- 3891 c 3892 c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2') 3893 c CALL writeg1d(ngrid,1,tsurf,'tsurf','K') 3894 c CALL writeg1d(ngrid,1,ps,'ps','Pa') 3895 c CALL writeg1d(ngrid,nlayer,zt,'T','K') 3896 c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1') 3897 c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1') 3898 c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1') 3899 call writediagsoil(ngrid,"soiltemp","Soil temperature","K", 3900 & 3,tsoil(:,:,iflat)) 3901 do islope=1,nslope 3902 write(str2(1:2),'(i2.2)') islope 3903 call writediagsoil(ngrid,"soiltemp_slope"//str2, 3904 & "Soil temperature","K", 3905 & 3,tsoil(:,:,islope)) 3906 ENDDO 3907 3908 ! THERMALS STUFF 1D 3909 if(calltherm) then 3910 3911 call WRITEDIAGFI(ngrid,'lmax_th', 3912 & 'hauteur du thermique','point', 3913 & 0,lmax_th_out) 3914 call WRITEDIAGFI(ngrid,'zmax_th', 3915 & 'hauteur du thermique','m', 3916 & 0,zmax_th) 3917 call WRITEDIAGFI(ngrid,'hfmax_th', 3918 & 'maximum TH heat flux','K.m/s', 3919 & 0,hfmax_th) 3920 call WRITEDIAGFI(ngrid,'wstar', 3921 & 'maximum TH vertical velocity','m/s', 3922 & 0,wstar) 3923 3924 end if ! of if (calltherm) 3925 3926 call WRITEDIAGFI(ngrid,'w','vertical velocity' 3927 & ,'m/s',1,pw) 3928 call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2) 3929 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0, 3930 & tsurf(:,iflat)) 3931 do islope=1,nslope 3932 write(str2(1:2),'(i2.2)') islope 3933 call WRITEDIAGFI(ngrid,"tsurf_slope"//str2, 3934 & "Surface temperature","K",0, 3935 & tsurf(:,islope)) 3936 ENDDO 3937 call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu) 3938 call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv) 3939 3940 call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay) 3941 call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev) 3942 call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho) 3788 3789 !PREVIOUSLY IN 1D ONLY 3943 3790 call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", 3944 & "K.s-1",1,dtrad) 3945 call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate', 3946 & 'w.m-2',1,zdtsw) 3947 call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate', 3948 & 'w.m-2',1,zdtlw) 3949 call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness" 3950 & ,"kg.m-2",0,qsurf(:,igcm_co2,iflat)) 3951 do islope=1,nslope 3952 write(str2(1:2),'(i2.2)') islope 3953 call WRITEDIAGFI(ngrid,"co2ice_slope"//str2, 3954 & "co2 ice thickness" 3955 & ,"kg.m-2",0,qsurf(:,igcm_co2,islope)) 3956 ENDDO 3957 3958 if (igcm_co2.ne.0) then 3959 call co2sat(ngrid*nlayer,zt,zqsatco2) 3960 do ig=1,ngrid 3961 do l=1,nlayer 3962 satuco2(ig,l) = zq(ig,l,igcm_co2)* 3963 & (mmean(ig,l)/44.01)*zplay(ig,l)/zqsatco2(ig,l) 3964 3965 c write(*,*) "In PHYSIQMOD, pt,zt,time ",pt(ig,l) 3966 c & ,zt(ig,l),ptime 3967 enddo 3968 enddo 3969 endif 3970 3971 call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps) 3972 call WRITEDIAGFI(ngrid,'temp','Temperature ', 3973 & 'K JA',1,zt) 3974 c call WRITEDIAGFI(ngrid,'temp2','Temperature ', 3975 c & 'K JA2',1,pt) 3976 3977 c CALL writeg1d(ngrid,1,tau,'tau','SI') 3978 do iq=1,nq 3979 c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 3980 call WRITEDIAGFI(ngrid,trim(noms(iq)), 3981 & trim(noms(iq)),'kg/kg',1,zq(1,1,iq)) 3982 end do 3983 if (doubleq) then 3984 call WRITEDIAGFI(ngrid,'rdust','rdust', 3985 & 'm',1,rdust) 3986 endif ! doubleq 1D 3791 & "K.s-1",3,dtrad) 3792 3987 3793 if (rdstorm) then 3988 3794 call writediagfi(1,'aerosol_dust','opacity of env. dust','' 3989 & , 1,aerosol(:,:,iaer_dust_doubleq))3795 & ,3,aerosol(:,:,iaer_dust_doubleq)) 3990 3796 call writediagfi(1,'aerosol_stormdust', 3991 3797 & 'opacity of storm dust','' 3992 & , 1,aerosol(:,:,iaer_stormdust_doubleq))3798 & ,3,aerosol(:,:,iaer_stormdust_doubleq)) 3993 3799 call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion', 3994 & 'kg.m-2.s-1', 0,zdqsdif(1,igcm_dust_mass,iflat))3800 & 'kg.m-2.s-1',2,zdqsdif(:,igcm_dust_mass,iflat)) 3995 3801 do islope=1,nslope 3996 3802 write(str2(1:2),'(i2.2)') islope 3997 3803 call WRITEDIAGFI(ngrid,'dqsdifdustq_slope'//str2, 3998 3804 & 'diffusion', 3999 & 'kg.m-2.s-1', 0,zdqsdif(1,igcm_dust_mass,islope))3805 & 'kg.m-2.s-1',2,zdqsdif(:,igcm_dust_mass,islope)) 4000 3806 ENDDO 4001 3807 call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion', 4002 & 'kg.m-2.s-1', 0,zdqsdif(1,igcm_stormdust_mass,iflat))3808 & 'kg.m-2.s-1',2,zdqsdif(:,igcm_stormdust_mass,iflat)) 4003 3809 do islope=1,nslope 4004 3810 write(str2(1:2),'(i2.2)') islope 4005 3811 call WRITEDIAGFI(ngrid,'dqsdifrdsq_slope'//str2, 4006 3812 & 'diffusion', 4007 & 'kg.m-2.s-1', 0,zdqsdif(1,igcm_stormdust_mass,islope))3813 & 'kg.m-2.s-1',2,zdqsdif(:,igcm_stormdust_mass,islope)) 4008 3814 ENDDO 4009 call WRITEDIAGFI(ngrid,'mstormdtot', 4010 & 'total mass of stormdust only', 4011 & 'kg.m-2',0,mstormdtot) 4012 call WRITEDIAGFI(ngrid,'mdusttot', 4013 & 'total mass of dust only', 4014 & 'kg.m-2',0,mdusttot) 4015 call WRITEDIAGFI(ngrid,'tau_pref_scenario', 4016 & 'Prescribed dust ref opt depth at 610 Pa', 4017 & 'NU',0,tau_pref_scenario) 4018 call WRITEDIAGFI(ngrid,'tau_pref_gcm', 4019 & 'Dust ref opt depth at 610 Pa in the GCM', 4020 & 'NU',0,tau_pref_gcm) 4021 call WRITEDIAGFI(ngrid,'rdsdqsdust', 4022 & 'deposited surface stormdust mass', 4023 & 'kg.m-2.s-1',0,rdsdqdustsurf) 4024 call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr', 4025 & 'kg/kg',1,rdsqdust) 4026 call WRITEDIAGFI(ngrid,"stormfract", 4027 & "fraction of the mesh,with stormdust", 4028 & "none",0,totstormfract) 4029 call WRITEDIAGFI(ngrid,'rdsqsurf', 4030 & 'stormdust at surface', 4031 & 'kg.m-2',0,qsurf(:,igcm_stormdust_mass,iflat)) 4032 do islope=1,nslope 4033 write(str2(1:2),'(i2.2)') islope 4034 call WRITEDIAGFI(ngrid,'rdsqsurf_slope'//str2, 4035 & 'stormdust at surface', 4036 & 'kg.m-2',0,qsurf(:,igcm_stormdust_mass,islope)) 4037 ENDDO 4038 call WRITEDIAGFI(ngrid,'qsurf', 4039 & 'dust mass at surface', 4040 & 'kg.m-2',0,qsurf(:,igcm_dust_mass,iflat)) 4041 do islope=1,nslope 4042 write(str2(1:2),'(i2.2)') islope 4043 call WRITEDIAGFI(ngrid,'qsurf_slope'//str2, 4044 & 'dust mass at surface', 4045 & 'kg.m-2',0,qsurf(:,igcm_dust_mass,islope)) 4046 ENDDO 4047 call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust', 4048 & 'm/s',1,wspeed) 4049 call WRITEDIAGFI(ngrid,'totaldustq','total dust mass', 4050 & 'kg/kg',1,qdusttotal) 4051 call WRITEDIAGFI(ngrid,'dsords', 4052 & 'density scaled opacity of stormdust', 4053 & 'm2.kg-1',1,dsords) 4054 call WRITEDIAGFI(ngrid,'zdqsed_dustq' 4055 & ,'sedimentation q','kg.m-2.s-1',1, 4056 & zdqsed(1,:,igcm_dust_mass)) 4057 call WRITEDIAGFI(ngrid,'zdqsed_rdsq' 4058 & ,'sedimentation q','kg.m-2.s-1',1, 4059 & zdqsed(1,:,igcm_stormdust_mass)) 4060 endif !(rdstorm 1D) 4061 4062 if (water.AND.tifeedback) then 4063 call WRITEDIAGFI(ngrid,"soiltemp", 4064 & "Soil temperature","K", 4065 & 1,tsoil) 4066 call WRITEDIAGFI(ngrid,'soilti', 4067 & 'Soil Thermal Inertia', 4068 & 'J.s-1/2.m-2.K-1',1,inertiesoil) 4069 endif 4070 4071 cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc 4072 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 4073 IF (water) THEN 4074 4075 if (.not.activice) then 4076 4077 tauTES=0 4078 do l=1,nlayer 4079 Qabsice = min( 4080 & max(0.4e6*rice(1,l)*(1.+nuice_ref)-0.05 ,0.),1.2 4081 & ) 4082 opTES(1,l)= 0.75 * Qabsice * 4083 & zq(1,l,igcm_h2o_ice) * 4084 & (zplev(1,l) - zplev(1,l+1)) / g 4085 & / (rho_ice * rice(1,l) * (1.+nuice_ref)) 4086 tauTES=tauTES+ opTES(1,l) 4087 enddo 4088 CALL WRITEDIAGFI(ngrid,'tauTESap', 4089 & 'tau abs 825 cm-1', 4090 & '',0,tauTES) 4091 else 4092 4093 CALL WRITEDIAGFI(ngrid,'tauTES', 4094 & 'tau abs 825 cm-1', 4095 & '',0,taucloudtes) 4096 endif 4097 4098 mtot = 0 4099 icetot = 0 4100 h2otot = qsurf(1,igcm_h2o_ice,iflat) 4101 if (hdo) THEN 4102 mtotD = 0 4103 icetotD = 0 4104 hdotot = qsurf(1,igcm_hdo_ice,iflat) 4105 ENDIF !hdo 4106 4107 do l=1,nlayer 4108 mtot = mtot + zq(1,l,igcm_h2o_vap) 4109 & * (zplev(1,l) - zplev(1,l+1)) / g 4110 icetot = icetot + zq(1,l,igcm_h2o_ice) 4111 & * (zplev(1,l) - zplev(1,l+1)) / g 4112 if (hdo) THEN 4113 mtotD = mtotD + zq(1,l,igcm_hdo_vap) 4114 & * (zplev(1,l) - zplev(1,l+1)) / g 4115 icetotD = icetotD + zq(1,l,igcm_hdo_ice) 4116 & * (zplev(1,l) - zplev(1,l+1)) / g 4117 ENDIF !hdo 4118 end do 4119 h2otot = h2otot+mtot+icetot 4120 IF (hdo) then 4121 hdotot = hdotot+mtotD+icetotD 4122 ENDIF ! hdo 4123 4124 4125 CALL WRITEDIAGFI(ngrid,'h2otot', 4126 & 'h2otot', 4127 & 'kg/m2',0,h2otot) 4128 CALL WRITEDIAGFI(ngrid,'mtot', 4129 & 'mtot', 4130 & 'kg/m2',0,mtot) 4131 CALL WRITEDIAGFI(ngrid,'icetot', 4132 & 'icetot', 4133 & 'kg/m2',0,icetot) 4134 CALL WRITEDIAGFI(ngrid,'albedo', 4135 & 'albedo', 4136 & '',2,albedo(1,1,iflat)) 4137 do islope=1,nslope 4138 write(str2(1:2),'(i2.2)') islope 4139 CALL WRITEDIAGFI(ngrid,'albedo_slope'//str2, 4140 & 'albedo', 4141 & '',2,albedo(1,1,islope)) 4142 ENDDO 4143 IF (hdo) THEN 4144 CALL WRITEDIAGFI(ngrid,'mtotD', 4145 & 'mtotD', 4146 & 'kg/m2',0,mtotD) 4147 CALL WRITEDIAGFI(ngrid,'icetotD', 4148 & 'icetotD', 4149 & 'kg/m2',0,icetotD) 4150 CALL WRITEDIAGFI(ngrid,'hdotot', 4151 & 'hdotot', 4152 & 'kg/m2',0,hdotot) 4153 4154 C Calculation of the D/H ratio 4155 do l=1,nlayer 4156 if (zq(1,l,igcm_h2o_vap).gt.qperemin) then 4157 DoH_vap(1,l) = 0.5*( zq(1,l,igcm_hdo_vap)/ 4158 & zq(1,l,igcm_h2o_vap) )/155.76e-6 4159 else 4160 DoH_vap(1,l) = 0. 4161 endif 4162 enddo 4163 4164 do l=1,nlayer 4165 if (zq(1,l,igcm_h2o_ice).gt.qperemin) then 4166 DoH_ice(1,l) = 0.5*( zq(1,l,igcm_hdo_ice)/ 4167 & zq(1,l,igcm_h2o_ice) )/155.76e-6 4168 else 4169 DoH_ice(1,l) = 0. 4170 endif 4171 enddo 4172 4173 CALL WRITEDIAGFI(ngrid,'DoH_vap', 4174 & 'D/H ratio in vapor', 4175 & ' ',1,DoH_vap) 4176 CALL WRITEDIAGFI(ngrid,'DoH_ice', 4177 & 'D/H ratio in ice', 4178 & '',1,DoH_ice) 4179 4180 ENDIF !Hdo 4181 4182 4183 if (scavenging) then 4184 4185 rave = 0 4186 do l=1,nlayer 4187 cccc Column integrated effective ice radius 4188 cccc is weighted by total ice surface area (BETTER) 4189 rave = rave + tauscaling(1) * 4190 & zq(1,l,igcm_ccn_number) * 4191 & (zplev(1,l) - zplev(1,l+1)) / g * 4192 & rice(1,l) * rice(1,l)* (1.+nuice_ref) 4193 enddo 4194 rave=icetot*0.75/max(rave*pi*rho_ice,1.e-30) ! surface weight 4195 4196 Nccntot= 0 4197 call watersat(ngrid*nlayer,zt,zplay,zqsat) 4198 do l=1,nlayer 4199 Nccntot = Nccntot + 4200 & zq(1,l,igcm_ccn_number)*tauscaling(1) 4201 & *(zplev(1,l) - zplev(1,l+1)) / g 4202 satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l) 4203 satu(1,l) = (max(satu(1,l),float(1))-1) 4204 ! & * zq(1,l,igcm_h2o_vap) * 4205 ! & (zplev(1,l) - zplev(1,l+1)) / g 4206 enddo 4207 call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1, 4208 & satu) 4209 CALL WRITEDIAGFI(ngrid,'Nccntot', 4210 & 'Nccntot', 4211 & 'nbr/m2',0,Nccntot) 4212 4213 call WRITEDIAGFI(ngrid,'zdqsed_dustq' 4214 & ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass)) 4215 call WRITEDIAGFI(ngrid,'zdqsed_dustN' 4216 &,'sedimentation N','Nbr.m-2.s-1',1, 4217 & zdqsed(1,:,igcm_dust_number)) 4218 4219 else ! of if (scavenging) 4220 4221 cccc Column integrated effective ice radius 4222 cccc is weighted by total ice mass (LESS GOOD) 4223 rave = 0 4224 do l=1,nlayer 4225 rave = rave + zq(1,l,igcm_h2o_ice) 4226 & * (zplev(1,l) - zplev(1,l+1)) / g 4227 & * rice(1,l) * (1.+nuice_ref) 4228 enddo 4229 rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight 4230 endif ! of if (scavenging) 4231 4232 4233 CALL WRITEDIAGFI(ngrid,'reffice', 4234 & 'reffice', 4235 & 'm',0,rave) 4236 4237 !Alternative A. Pottier weighting 4238 rave2 = 0. 4239 totrave2 = 0. 4240 do l=1,nlayer 4241 rave2 =rave2+ zq(1,l,igcm_h2o_ice)*rice(1,l) 4242 totrave2 = totrave2 + zq(1,l,igcm_h2o_ice) 4243 end do 4244 rave2=max(rave2/max(totrave2,1.e-30),1.e-30) 4245 CALL WRITEDIAGFI(ngrid,'rmoym', 4246 & 'reffice', 4247 & 'm',0,rave2) 4248 4249 do iq=1,nq 4250 call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s', 4251 & trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq,iflat)) 4252 end do 4253 4254 call WRITEDIAGFI(ngrid,"watercap","Water ice thickness" 4255 & ,"kg.m-2",0,watercap) 3815 endif !(rdstorm) 3816 3817 if(water) then 3818 if (.not.scavenging) then 4256 3819 call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice', 4257 3820 & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)) … … 4266 3829 call WRITEDIAGFI(ngrid,'zdqcloud_vapD','cloud vap hdo', 4267 3830 & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_vap)) 4268 4269 3831 ENDIF ! hdo 4270 4271 call WRITEDIAGFI(ngrid,"rice","ice radius","m",1, 4272 & rice) 4273 4274 if (CLFvarying) then 4275 call WRITEDIAGFI(ngrid,'totcloudfrac', 4276 & 'Total cloud fraction', 4277 & ' ',0,totcloudfrac) 4278 endif !clfvarying 4279 3832 endif !not.scavenging 4280 3833 ENDIF ! of IF (water) 4281 4282 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 4283 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 4284 4285 zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g 4286 4287 do l=2,nlayer-1 4288 tmean=zt(1,l) 4289 if(zt(1,l).ne.zt(1,l-1)) 4290 & tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1)) 4291 zlocal(l)= zlocal(l-1) 4292 & -log(zplay(1,l)/zplay(1,l-1))*rnew(1,l)*tmean/g 4293 enddo 4294 zlocal(nlayer)= zlocal(nlayer-1)- 4295 & log(zplay(1,nlayer)/zplay(1,nlayer-1))* 4296 & rnew(1,nlayer)*tmean/g 4297 4298 3834 !PREVIOUSLY IN 1D ONLY 3835 3836 c ========================================================== 3837 c END OF WRITEDIAGFI 3838 c ========================================================== 4299 3839 #endif 4300 4301 END IF ! if(ngrid.ne.1) 3840 ! of ifdef MESOSCALE 3841 3842 c ELSE ! if(ngrid.eq.1) 3843 3844 c#ifndef MESOSCALE 3845 c write(*, 3846 c & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)') 3847 c & zls*180./pi,odpref,tau_pref_scenario 3848 c#endif 3849 3850 c END IF ! if(ngrid.ne.1) 3851 3852 4302 3853 ! test for co2 conservation with co2 microphysics 4303 3854 if (igcm_co2_ice.ne.0) then -
trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90
r2930 r2931 642 642 ! WRITEDIAGFI 643 643 ! ********************************************************************** 644 IF (ngrid.eq.1) THEN645 CALL WRITEDIAGFI(ngrid,'wup_top', &646 'wup_top','',1,wup(:,:))647 CALL WRITEDIAGFI(ngrid,'wrad_top', &648 'wrad_top','',1,wrad(:,:))649 CALL WRITEDIAGFI(ngrid,'wfin_top', &650 'wfin_top','',1,wfin(:,:))651 CALL WRITEDIAGFI(ngrid,'wlwmax_top', &652 'wlwmax_top','',0,wup(:,lwmax(1)))653 CALL WRITEDIAGFI(ngrid,'w0_top', &654 'w0_top','',0,wup(:,lsummit(1)+1))655 CALL WRITEDIAGFI(ngrid,'alpha_hmons', &656 'alpha_hmons','',0,alpha_hmons)657 CALL WRITEDIAGFI(ngrid,'zt_top', &658 'zt_top','',1,t_top(:,:))659 CALL WRITEDIAGFI(ngrid,'dt_top', &660 'dt_top','',0,dt_top(:))661 CALL WRITEDIAGFI(ngrid,'envtemp', &662 'envtemp','',1,zt(:,:))663 CALL WRITEDIAGFI(ngrid,'t_env', &664 't_env','',1,t_env(:,:))665 CALL WRITEDIAGFI(ngrid,'newzt_top', &666 'newzt_top','',1,newzt(:,:))667 CALL WRITEDIAGFI(ngrid,'deltahr_top', &668 'deltahr_top','',1,deltahr(:,:))669 CALL WRITEDIAGFI(ngrid,'rholev', &670 'rholev','',1,rho(:,:))671 CALL WRITEDIAGFI(ngrid,'rhobarz', &672 'rhobarz','',1,rhobarz(:,:))673 CALL WRITEDIAGFI(ngrid,'zlsummit', &674 'zlsummit','',0,pzlay(:,lsummit(1)))675 CALL WRITEDIAGFI(ngrid,'zlwmax', &676 'zlwmax','',0,pzlay(:,lwmax(1)))677 CALL WRITEDIAGFI(ngrid,'pzlev_top', &678 'pzlev_top','',1,pzlev(:,:))679 CALL WRITEDIAGFI(ngrid,'coefdetrain', &680 'coefdetrain','',1,coefdetrain(:,:))681 CALL WRITEDIAGFI(ngrid,'zdtvert', &682 'zdtvert','',1,zdtvert(:,:))683 CALL WRITEDIAGFI(ngrid,'entr', &684 'entr','',0,entr(:))685 ELSE686 644 CALL WRITEDIAGFI(ngrid,'wup_top', & 687 645 'wup_top','',3,wup(:,:)) 646 CALL WRITEDIAGFI(ngrid,'wlwmax_top', & 647 'wlwmax_top','',2,wup(:,lwmax(1))) 648 CALL WRITEDIAGFI(ngrid,'w0_top', & 649 'w0_top','',2,wup(:,lsummit(1)+1)) 688 650 CALL WRITEDIAGFI(ngrid,'wup2_top', & 689 651 'wup2_top','',3,wup2(:,:)) … … 700 662 CALL WRITEDIAGFI(ngrid,'zt_top', & 701 663 'zt_top','',3,t_top(:,:)) 664 CALL WRITEDIAGFI(ngrid,'dt_top', & 665 'dt_top','',2,dt_top(:)) 702 666 CALL WRITEDIAGFI(ngrid,'ztemp', & 703 667 'envtemp','',3,zt(:,:)) … … 719 683 'detrain_rate', & 720 684 '',3,detrain_rate(:,:)) 721 722 ENDIF 685 CALL WRITEDIAGFI(ngrid,'rholev', & 686 'rholev','',3,rho(:,:)) 687 CALL WRITEDIAGFI(ngrid,'zlsummit', & 688 'zlsummit','',2,pzlay(:,lsummit(1))) 689 CALL WRITEDIAGFI(ngrid,'pzlev_top', & 690 'pzlev_top','',3,pzlev(:,:)) 691 CALL WRITEDIAGFI(ngrid,'entr', & 692 'entr','',2,entr(:)) 723 693 724 694 END SUBROUTINE topmons
Note: See TracChangeset
for help on using the changeset viewer.
