Changeset 2931 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Apr 3, 2023, 4:49:33 PM (20 months ago)
Author:
romain.vande
Message:

Mars PCM:
Remove the test ngrid==1 for writediagfi. The routine can handle well the 1D case.
RV

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r2926 r2931  
    25122512
    25132513#ifndef MESOSCALE
    2514             IF (ngrid .eq. 1) THEN
    2515                dimout=0
    2516             ELSE
    2517                dimout=2
    2518             ENDIF
    25192514            DO n=1,n_out
    25202515               write(zstring, '(F8.6)') z_out(n)
    25212516               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))
    25232518               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))
    25252520            ENDDO
    25262521            call WRITEDIAGFI(ngrid,'u_star',
    2527      &   'friction velocity','m/s',dimout,ustar)
     2522     &   'friction velocity','m/s',2,ustar)
    25282523            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)
    25322525            call WRITEDIAGFI(ngrid,'vvv',
    2533      &   'Vertical velocity variance at zout','m',dimout,vvv)
     2526     &   'Vertical velocity variance at zout','m',2,vvv)
    25342527            call WRITEDIAGFI(ngrid,'vhf',
    2535      &   'Vertical heat flux at zout','m',dimout,vhf)
     2528     &   'Vertical heat flux at zout','m',2,vhf)
    25362529#else
    25372530         T_out1(:)=T_out(:,1)
    … …  
    26102603#endif
    26112604
    2612          IF (ngrid.NE.1) then
     2605c         IF (ngrid.NE.1) then
    26132606
    26142607c        -------------------------------------------------------------------
    … …  
    26402633           if(doubleq) then
    26412634              do ig=1,ngrid
     2635         IF (sedimentation) THEN
    26422636                dqdustsurf(ig) =
    26432637     &                zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
    26442638                dndustsurf(ig) =
    26452639     &                zdqssed(ig,igcm_dust_number)*tauscaling(ig)
     2640         ENDIF
    26462641                ndust(ig,:) =
    26472642     &                zq(ig,:,igcm_dust_number)*tauscaling(ig)
    … …  
    26512646              if (scavenging) then
    26522647                do ig=1,ngrid
     2648         IF (sedimentation) THEN
    26532649                  dqdustsurf(ig) = dqdustsurf(ig) +
    26542650     &                     zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
    26552651                  dndustsurf(ig) = dndustsurf(ig) +
    26562652     &                     zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
     2653         ENDIF
    26572654                  nccn(ig,:) =
    26582655     &                     zq(ig,:,igcm_ccn_number)*tauscaling(ig)
    … …  
    28532850          call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
    28542851     &                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)
    28612852          call wstats(ngrid,"fluxsurf_dir_dn_sw",
    28622853     &                "Direct incoming SW flux at surface",
    … …  
    29502941             
    29512942             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))
    29622943               call wstats(ngrid,'dqsdust',
    29632944     &                        'deposited surface dust mass',
    … …  
    31593140     &    "Surface emissivity","w.m-1",2, emis(:,islope))
    31603141         ENDDO
    3161 c        call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
    3162 c        call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
    31633142         call WRITEDIAGFI(ngrid,"zzlev","Interlayer altitude",
    31643143     &                    "m",3,zzlev(:,1:nlayer))
    … …  
    32233202         call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
    32243203         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)
    32273204         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)
    32303206        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
    32313207     &                   'K/s',3,zdtsw)
    … …  
    35013477        end if !(water)
    35023478
    3503         if (water.and..not.photochem) then
    3504           iq = nq
    3505 c          write(str2(1:2),'(i2.2)') iq
    3506 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 
    35183479c        ----------------------------------------------------------
    35193480c        Outputs of the dust cycle
    … …  
    35403501#ifndef MESOINI
    35413502           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))
    35563503             call WRITEDIAGFI(ngrid,'dqsdust',
    35573504     &                        'deposited surface dust mass',
    … …  
    36483595     &          ,'sedimentation q','kg.m-2.s-1',3,
    36493596     &          zdqsed(:,:,igcm_stormdust_mass))
     3597             call WRITEDIAGFI(ngrid,'zdqsed_dustN'
     3598     &          ,'sedimentation N','Nbr.m-2.s-1',3,
     3599     &          zdqsed(:,:,igcm_dust_number))
    36503600             call WRITEDIAGFI(ngrid,'rdust','rdust',
    36513601     &                        'm',3,rdust)
    … …  
    37143664         ENDDO
    37153665           endif ! (scavenging)
    3716 
    3717 c          if (submicron) then
    3718 c            call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
    3719 c    &                        'kg/kg',3,pq(1,1,igcm_dust_submicron))
    3720 c          endif ! (submicron)
    37213666
    37223667#else
    … …  
    37893734     $           0,[PhiEscD])
    37903735
    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 
    38073736         endif  !(callthermos)
    38083737
    … …  
    38223751c        ----------------------------------------------------------
    38233752      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)
    38393756        call WRITEDIAGFI(ngrid,'zmax_th',
    38403757     &                   'hauteur du thermique','m',
    … …  
    38693786     &                     3,tsoil(:,:,islope))
    38703787         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
    39433790         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
    39873793           if (rdstorm) then
    39883794             call writediagfi(1,'aerosol_dust','opacity of env. dust',''
    3989      &                        ,1,aerosol(:,:,iaer_dust_doubleq))
     3795     &                        ,3,aerosol(:,:,iaer_dust_doubleq))
    39903796             call writediagfi(1,'aerosol_stormdust',
    39913797     &                         'opacity of storm dust',''
    3992      &                        ,1,aerosol(:,:,iaer_stormdust_doubleq))
     3798     &                        ,3,aerosol(:,:,iaer_stormdust_doubleq))
    39933799             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))
    39953801         do islope=1,nslope
    39963802          write(str2(1:2),'(i2.2)') islope
    39973803             call WRITEDIAGFI(ngrid,'dqsdifdustq_slope'//str2,
    39983804     &             '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))
    40003806         ENDDO
    40013807             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))
    40033809         do islope=1,nslope
    40043810          write(str2(1:2),'(i2.2)') islope
    40053811             call WRITEDIAGFI(ngrid,'dqsdifrdsq_slope'//str2,
    40063812     &           '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))
    40083814         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
    42563819        call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
    42573820     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
    … …  
    42663829        call WRITEDIAGFI(ngrid,'zdqcloud_vapD','cloud vap hdo',
    42673830     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_vap))
    4268 
    42693831        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
    42803833        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
     3836c        ==========================================================
     3837c        END OF WRITEDIAGFI
     3838c        ==========================================================
    42993839#endif
    4300 
    4301       END IF       ! if(ngrid.ne.1)
     3840! of ifdef MESOSCALE
     3841
     3842c      ELSE     ! if(ngrid.eq.1)
     3843
     3844c#ifndef MESOSCALE
     3845c         write(*,
     3846c     &    '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
     3847c     &    zls*180./pi,odpref,tau_pref_scenario
     3848c#endif
     3849
     3850c      END IF       ! if(ngrid.ne.1)
     3851
     3852
    43023853! test for co2 conservation with co2 microphysics
    43033854      if (igcm_co2_ice.ne.0) then

  • trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90

    r2930 r2931  
    642642! WRITEDIAGFI
    643643! **********************************************************************
    644         IF (ngrid.eq.1) THEN
    645                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         ELSE
    686644               CALL WRITEDIAGFI(ngrid,'wup_top', &
    687645                '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))
    688650               CALL WRITEDIAGFI(ngrid,'wup2_top', &
    689651                'wup2_top','',3,wup2(:,:))
    … …  
    700662               CALL WRITEDIAGFI(ngrid,'zt_top', &
    701663                'zt_top','',3,t_top(:,:))
     664               CALL WRITEDIAGFI(ngrid,'dt_top', &
     665                'dt_top','',2,dt_top(:))
    702666               CALL WRITEDIAGFI(ngrid,'ztemp', &
    703667                'envtemp','',3,zt(:,:))
    … …  
    719683                              'detrain_rate', &
    720684                              '',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(:))
    723693       
    724694        END SUBROUTINE topmons
Note: See TracChangeset for help on using the changeset viewer.