Ignore:
Timestamp:
Apr 4, 2018, 4:18:35 PM (7 years ago)
Author:
jvatant
Message:

Ooops just realized that in r1908 physiq_mod.F90
hasn't been pushed ...
+ Also add "cosmetic" changes for arrays (1:nlayer)->(:)
--JVO

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

    r1904 r1914  
    3636      use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
    3737      use time_phylmdz_mod, only: daysec
    38       use logic_mod, only: moyzon_ch, moyzon_mu
     38      use logic_mod, only: moyzon_ch
    3939      use moyzon_mod, only: tmoy, playmoy, zphibar, zphisbar, zplevbar, &
    4040                            zplaybar, zzlevbar, zzlaybar, ztfibar, zqfibar
    … …  
    367367! ----------------------------------------------------
    368368 
    369       real ctimestep ! Chemistry timestep (s)
     369      real,save :: ctimestep ! Chemistry timestep (s)
     370!$OMP THREADPRIVATE(ctimestep)
    370371 
    371372      ! Chemical tracers in molar fraction
    372       real, dimension(:,:,:), allocatable, save   :: ychim ! (mol/mol)
    373 !$OMP THREADPRIVATE(ychim)
     373      real, dimension(ngrid,nlayer,nkim)          :: ychim ! (mol/mol)
    374374
    375375      ! Molar fraction tendencies ( chemistry and condensation ) for tracers (mol/mol/s)
    … …  
    381381      real, dimension(:,:), allocatable, save     :: qysat ! (mol/mol)
    382382!$OMP THREADPRIVATE(qysat)
    383       real temp_eq(nlayer), press_eq(nlayer) ! Planetary averages for the init. of saturation profiles.
     383      real temp_eq(nlayer), press_eq(nlayer) ! Planetary averages for the init. of saturation profiles (K,mbar)
    384384
    385385      ! Surface methane tank
    … …  
    403403    !   Or one can put calmufi in MMP_GCM module (in muphytitan).
    404404    INTERFACE
    405       SUBROUTINE calmufi(dt, plev, zlev, play, zlay, temp, pq, zdq)
     405      SUBROUTINE calmufi(dt, plev, zlev, play, zlay, temp, pq, zdqfi, zdq)
    406406        REAL(kind=8), INTENT(IN)                 :: dt    !! Physics timestep (s).
    407407        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev  !! Pressure levels (Pa).
    … …  
    410410        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay  !! Altitude at the center of each layer (m).
    411411        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K). 
    412         REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq   !! Tracers (\(kg.kg^{-1}}\)).
    413         REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq  !! Tendency for tracers (\(kg.kg^{-1}}\)).
     412        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(kg.kg^{-1}}\)).
     413        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi !! Tendency from former processes for tracers (\(kg.kg^{-1}}\)).
     414        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq   !! Microphysical tendency for tracers (\(kg.kg^{-1}}\)).
    414415      END SUBROUTINE calmufi
    415416    END INTERFACE
    … …  
    505506         if ( callchim ) then
    506507
    507             allocate(ychim(ngrid,nlayer,nkim))
     508            if ( moyzon_ch .and. ngrid.eq.1 ) then
     509               print *, "moyzon_ch=",moyzon_ch," and ngrid=1"
     510               print *, "Please desactivate zonal mean for 1D !"
     511               stop
     512            endif
     513
    508514            allocate(dycchi(ngrid,nlayer,nkim)) ! only for chemical tracers
    509515            allocate(qysat(nlayer,nkim))
    … …  
    513519
    514520            ! qysat is taken at the equator ( small variations of t,p )
    515             temp_eq(:)  = tmoy(:)
    516             press_eq(:) = playmoy(:)/100. ! in mbar
     521            if (ngrid.ne.1) then ! TODO : a patcher (lecture d'un profil?) si jamais on a plus acces a moyzon_mod !
     522              temp_eq(:)  = tmoy(:)
     523              press_eq(:) = playmoy(:)/100. ! in mbar
     524            else
     525              temp_eq(:)  = pt(1,:)
     526              press_eq(:) = pplay(1,:)/100.
     527            endif
    517528           
    518529            call inicondens(press_eq,temp_eq,qysat)
    … …  
    553564           if (is_master) write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!"
    554565           do isoil=1,nsoilmx
    555              tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
     566             tsoil(:,isoil)=tsurf(:)
    556567           enddo
    557568           if (is_master) write(*,*) "Physiq: initializing day_ini to pdat !"
    … …  
    648659      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    649660
    650       pdt(1:ngrid,1:nlayer) = 0.0     
    651       zdtsurf(1:ngrid)      = 0.0
    652       pdq(1:ngrid,1:nlayer,1:nq) = 0.0
    653       dqsurf(1:ngrid,1:nq)= 0.0
    654       pdu(1:ngrid,1:nlayer) = 0.0
    655       pdv(1:ngrid,1:nlayer) = 0.0
    656       pdpsrf(1:ngrid)       = 0.0
    657       zflubid(1:ngrid)      = 0.0     
    658       taux(1:ngrid) = 0.0
    659       tauy(1:ngrid) = 0.0
     661      pdt(:,:)    = 0.0     
     662      zdtsurf(:)  = 0.0
     663      pdq(:,:,:) = 0.0
     664      dqsurf(:,:) = 0.0
     665      pdu(:,:)    = 0.0
     666      pdv(:,:)    = 0.0
     667      pdpsrf(:)   = 0.0
     668      zflubid(:)  = 0.0     
     669      taux(:)    = 0.0
     670      tauy(:)    = 0.0
    660671
    661672      zday=pday+ptime ! Compute time, in sols (and fraction thereof).
    … …  
    696707      ! Compute geopotential between layers.
    697708      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    698       zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
     709      zzlay(:,:)=pphi(:,:)
    699710      do l=1,nlayer         
    700          zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
     711         zzlay(:,l)= zzlay(:,l)/glat(:)
    701712      enddo
    702713
    703       zzlev(1:ngrid,1)=0.
    704       zzlev(1:ngrid,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
     714      zzlev(:,1)=0.
     715      zzlev(:,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
    705716
    706717      do l=2,nlayer
    … …  
    712723      enddo     
    713724
    714       ! Zonal averages needed for chemistry and microphysics
    715       ! ~~~~~~~~~~~~~ Taken from old Titan ~~~~~~~~~~~~~~~~~
    716       if (moyzon_ch .or. moyzon_mu) then
    717          
    718          print *, "------------------------------ </CRITICAL ALERT> -------------------------------"
    719          print *, "WARNING : YOU ARE USING ZONAL MEANS THAT SEEM TO LEAD TO CRASHES (-infinity) ..."
    720          print *, "------------------------------ <CRITICAL ALERT/> -------------------------------"
    721 
     725      ! Zonal averages needed for chemistry
     726      ! ~~~~~~~ Taken from old Titan ~~~~~~
     727      if (moyzon_ch) then
     728         
    722729         zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
    723730         ! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
    … …  
    770777     ! But first linearly interpolate mass flux to mid-layers
    771778      do l=1,nlayer-1
    772          pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
     779         pw(:,l)=0.5*(flxw(:,l)+flxw(:,l+1))
    773780      enddo
    774       pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
     781      pw(:,nlayer)=0.5*flxw(:,nlayer) ! since flxw(nlayer+1)=0
    775782      do l=1,nlayer
    776          pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /  &
    777                        (pplay(1:ngrid,l)*cell_area(1:ngrid))
     783         pw(:,l)=(pw(:,l)*r*pt(:,l)) / (pplay(:,l)*cell_area(:))
    778784      enddo
    779785
    … …  
    828834               ! computations... waste of time...
    829835               DO i = 1, ngrid
    830                  i2e(:) = ( pplev(i,1:nlayer)-pplev(i,2:nlayer+1) ) / g
     836                 i2e(1:nlayer) = ( pplev(i,1:nlayer)-pplev(i,2:nlayer+1) ) / g
    831837                 pqo(i,:,:) = 0.0
    832838                 DO j=1,nmicro-nice
    … …  
    847853               ! Radiative flux from the sky absorbed by the surface (W.m-2).
    848854               GSR=0.0
    849                fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid)
     855               fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)+fluxsurfabs_sw(:)
    850856
    851857                            !if(noradsurf)then ! no lower surface; SW flux just disappears
    852                             !   GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea
    853                             !   fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
     858                            !   GSR = SUM(fluxsurf_sw(:)*cell_area(:))/totarea
     859                            !   fluxrad_sky(:)=emis(:)*fluxsurf_lw(:)
    854860                            !   print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
    855861                            !endif
    856862
    857863               ! Net atmospheric radiative heating rate (K.s-1)
    858                dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
     864               dtrad(:,:)=zdtsw(:,:)+zdtlw(:,:)
    859865               
    860866            elseif(newtonian)then
    … …  
    865871               call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
    866872
    867                zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
     873               zdtsurf(:) = +(pt(:,1)-tsurf(:))/ptimestep
    868874               ! e.g. surface becomes proxy for 1st atmospheric layer ?
    869875
    … …  
    873879! II.c Atmosphere has no radiative effect
    874880! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    875                fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
     881               fluxtop_dn(:)  = fract(:)*mu0(:)*Fat1AU/dist_star**2
    876882               if(ngrid.eq.1)then ! / by 4 globally in 1D case...
    877883                  fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
    878884               endif
    879                fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
     885               fluxsurf_sw(:) = fluxtop_dn(:)
    880886               print*,'------------WARNING---WARNING------------' ! by MT2015.
    881887               print*,'You are in corrk=false mode, '
    882888               print*,'and the surface albedo is taken equal to the first visible spectral value'               
    883889               
    884                fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
    885                fluxrad_sky(1:ngrid)    = fluxsurfabs_sw(1:ngrid)
    886                fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
    887 
    888                dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
     890               fluxsurfabs_sw(:) = fluxtop_dn(:)*(1.-albedo(:,1))
     891               fluxrad_sky(:)    = fluxsurfabs_sw(:)
     892               fluxtop_lw(:)     = emis(:)*sigma*tsurf(:)**4
     893
     894               dtrad(:,:)=0.0 ! no atmospheric radiative heating
    889895
    890896            endif ! end of corrk
    … …  
    895901         ! Transformation of the radiative tendencies
    896902         ! ------------------------------------------
    897          zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
    898          zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid)
    899          fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
    900          pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
     903         zplanck(:)=tsurf(:)*tsurf(:)
     904         zplanck(:)=emis(:)*sigma*zplanck(:)*zplanck(:)
     905         fluxrad(:)=fluxrad_sky(:)-zplanck(:)
     906         pdt(:,:)=pdt(:,:)+dtrad(:,:)
    901907         
    902908         ! Test of energy conservation
    … …  
    931937      if (calldifv) then
    932938     
    933          zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
     939         zflubid(:)=fluxrad(:)+fluxgrd(:)
    934940
    935941         ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
    … …  
    947953         else
    948954         
    949             zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
     955            zdh(:,:)=pdt(:,:)/zpopsk(:,:)
    950956 
    951957            call vdifc(ngrid,nlayer,nq,zpopsk,                &
    … …  
    958964                       taux,tauy,lastcall)
    959965
    960             zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
    961             zdqevap(1:ngrid,1:nlayer)=0.
     966            zdtdif(:,:)=zdhdif(:,:)*zpopsk(:,:) ! for diagnostic only
     967            zdqevap(:,:)=0.
    962968
    963969         end if !end of 'UseTurbDiff'
    964970
    965971
    966          pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
    967          pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
    968          pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
    969          zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
     972         pdv(:,:)=pdv(:,:)+zdvdif(:,:)
     973         pdu(:,:)=pdu(:,:)+zdudif(:,:)
     974         pdt(:,:)=pdt(:,:)+zdtdif(:,:)
     975         zdtsurf(:)=zdtsurf(:)+zdtsdif(:)
    970976
    971977         if (tracer) then
    972            pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
    973            dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
     978           pdq(:,:,:)=pdq(:,:,:)+ zdqdif(:,:,:)
     979           dqsurf(:,:)=dqsurf(:,:) + zdqsdif(:,:)
    974980         end if ! of if (tracer)
    975981
    … …  
    10101016         if(.not.newtonian)then
    10111017
    1012             zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
     1018            zdtsurf(:) = zdtsurf(:) + (fluxrad(:) + fluxgrd(:))/capcal(:)
    10131019
    10141020         endif
    … …  
    10231029      if(calladj) then
    10241030
    1025          zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
    1026          zduadj(1:ngrid,1:nlayer)=0.0
    1027          zdvadj(1:ngrid,1:nlayer)=0.0
    1028          zdhadj(1:ngrid,1:nlayer)=0.0
     1031         zdh(:,:) = pdt(:,:)/zpopsk(:,:)
     1032         zduadj(:,:)=0.0
     1033         zdvadj(:,:)=0.0
     1034         zdhadj(:,:)=0.0
    10291035
    10301036
    … …  
    10361042                      zdqadj)
    10371043
    1038          pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
    1039          pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer)
    1040          pdt(1:ngrid,1:nlayer)    = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
    1041          zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
     1044         pdu(:,:) = pdu(:,:) + zduadj(:,:)
     1045         pdv(:,:) = pdv(:,:) + zdvadj(:,:)
     1046         pdt(:,:)    = pdt(:,:) + zdhadj(:,:)*zpopsk(:,:)
     1047         zdtadj(:,:) = zdhadj(:,:)*zpopsk(:,:) ! for diagnostic only
    10421048
    10431049         if(tracer) then
    1044             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
     1050            pdq(:,:,:) = pdq(:,:,:) + zdqadj(:,:,:)
    10451051         end if
    10461052
    … …  
    10601066
    10611067      if (tracer) then
    1062      
     1068
    10631069  ! -------------------------
    10641070  !   V.1. Chemistry
    … …  
    10671073         if (callchim) then
    10681074
    1069             ! Convert to molar fraction
    1070             if (moyzon_ch) then
    1071               do iq = 1,nkim
    1072                  ychim(:,:,iq) = zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
    1073               enddo
    1074             else
    1075               do iq = 1,nkim
    1076                  ychim(:,:,iq) = pq(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
    1077               enddo
    1078             endif
    1079 
    1080             ! Condensation tendency after the transport
     1075            ! o. Convert updated tracers to molar fraction
     1076            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     1077            do iq = 1,nkim
     1078               ychim(:,:,iq) = ( pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro) ) / rat_mmol(iq+nmicro)
     1079            enddo
     1080
     1081            ! i. Condensation after the transport
     1082            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    10811083            do iq=1,nkim
    10821084               do l=1,nlayer
    10831085                  do ig=1,ngrid
    10841086                     if ( ychim(ig,l,iq).gt.qysat(l,iq) ) then
    1085                         dyccond(ig,l,iq+nmicro)= ( -ychim(ig,l,iq)+qysat(l,iq) ) / ptimestep
     1087                        dyccond(ig,l,iq+nmicro) = ( -ychim(ig,l,iq)+qysat(l,iq) ) / ptimestep
     1088                     else
     1089                        dyccond(ig,l,iq+nmicro) = 0.0 ! since array not saved
    10861090                     endif
    10871091                  enddo
    … …  
    10891093            enddo
    10901094
    1091             if ( callclouds ) dyccond(:,:,ices_indx) = 0.0 ! Condensation will be calculated in the cloud microphysics
    1092 
     1095            if ( callclouds ) dyccond(:,:,ices_indx) = 0.0 ! Condensation will be calculated in the cloud microphysics
     1096
     1097            do iq=1,nkim
     1098              ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro) ! update molar ychim for following calchim
     1099              zdqcond(:,:,iq+nmicro) = dyccond(:,:,iq+nmicro)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
     1100            enddo
     1101           
     1102            pdq(:,:,:) = pdq(:,:,:) + zdqcond(:,:,:)
     1103
     1104            ! ii. Photochemistry
     1105            ! ~~~~~~~~~~~~~~~~~~
    10931106            if( mod(icount-1,ichim).eq.0. ) then
    10941107               
    10951108               print *, "We enter in the chemistry ..."
    10961109
    1097                if (moyzon_ch.and.(ngrid.ne.1)) then ! 2D zonally averaged chemistry !! SEEMS TO CRASH DON'T USE IT !!
    1098                  call calchim(ngrid,ychim,declin,zls,ctimestep,ztfibar,zphibar,  &
     1110               if (moyzon_ch) then ! 2D zonally averaged chemistry
     1111
     1112                 do iq = 1,nkim ! In this case we send zonal average from dynamics to chem. module
     1113                   ychim(:,:,iq) = zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
     1114                 enddo
     1115
     1116                 call calchim(ngrid,ychim,declin,ctimestep,ztfibar,zphibar,  &
    10991117                              zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
     1118
    11001119               else ! 3D chemistry (or 1D run)
    1101                  call calchim(ngrid,ychim,declin,zls,ctimestep,pt,pphi,  &
     1120                 call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,  &
    11021121                              pplay,pplev,zzlay,zzlev,dycchi)
    1103                endif
    1104 
     1122               endif ! if moyzon
    11051123            endif
    11061124           
    1107             ! TEMPORARY COMMENT
    1108             ! These conversions as well as 2D/1D should be put in phytrac
    1109             ! ( GCM-friendly tracers and tendencies -> format for photochem and microphys )
    1110 
    1111             ! We convert tendencies to mass mixing ratio
    11121125            do iq=1,nkim
    1113                zdqchi(:,:,iq+nmicro)  = dycchi(:,:,iq) * rat_mmol(iq+nmicro)
    1114                zdqcond(:,:,iq+nmicro) = dyccond(:,:,iq+nmicro) * rat_mmol(iq+nmicro)
    1115             enddo
    1116 
    1117             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
    1118                  zdqchi(1:ngrid,1:nlayer,1:nq) + zdqcond(1:ngrid,1:nlayer,1:nq)
    1119                
     1126               zdqchi(:,:,iq+nmicro) = dycchi(:,:,iq)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
     1127 
     1128               where( (pq(:,:,iq+nmicro)+zdqchi(:,:,iq+nmicro) ).LT.1e-40)  & ! When using zonal means we set the same tendency
     1129                       zdqchi(:,:,iq+nmicro) = 1e-40 - ( pq(:,:,iq+nmicro) )  ! everywhere in longitude -> can lead to negs !
     1130             enddo
     1131
     1132            pdq(:,:,:) = pdq(:,:,:) + zdqchi(:,:,:)
    11201133           
    11211134         endif ! end of 'callchim'
    … …  
    11271140         if (callmufi) then
    11281141#ifdef USE_QTEST
    1129                call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,pt,tpq,zdqmufi)
    1130                tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimesep here
     1142            dtpq(:,:,:) = 0.0 ! we want tpq to go only through mufi
     1143            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,pt,tpq,dtpq,zdqmufi)
     1144            tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
    11311145#else
    1132             ! Inside this routine we will split 2D->1D, intensive->extensive and separate different types of tracers
    1133             ! Should be put in phytrac
    1134 
    1135                if (moyzon_mu.and.(ngrid.ne.1)) then ! 2D zonally averaged microphysics !! SEEMS TO CRASH DON'T USE IT !!
    1136                   call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar,ztfibar,zqfibar,zdqmufi)
    1137                else ! 3D microphysics (or 1D run)
    1138                   call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,pt,pq,zdqmufi)
    1139                endif
    1140 
    1141             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmufi(1:ngrid,1:nlayer,1:nq)
     1146            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,pt,pq,pdq,zdqmufi)
     1147            pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
    11421148#endif
    1143          endif ! end of 'callmufi'
     1149         endif
    11441150     
    11451151  ! ---------------
    … …  
    11501156         if(mass_redistrib) then
    11511157
    1152             zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * zdqevap(1:ngrid,1:nlayer)
     1158            zdmassmr(:,:) = mass(:,:) * zdqevap(:,:)
    11531159
    11541160            do ig = 1, ngrid
    1155                zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
     1161               zdmassmr_col(ig)=SUM(zdmassmr(ig,:))
    11561162            enddo
    11571163           
    … …  
    11651171                                     zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
    11661172         
    1167             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
    1168             dqsurf(1:ngrid,1:nq)       = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
    1169             pdt(1:ngrid,1:nlayer)      = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
    1170             pdu(1:ngrid,1:nlayer)      = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
    1171             pdv(1:ngrid,1:nlayer)      = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
    1172             pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
    1173             zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
     1173            pdq(:,:,:)  = pdq(:,:,:)  + zdqmr(:,:,:)
     1174            dqsurf(:,:) = dqsurf(:,:) + zdqsurfmr(:,:)
     1175            pdt(:,:)    = pdt(:,:)    + zdtmr(:,:)
     1176            pdu(:,:)    = pdu(:,:)    + zdumr(:,:)
     1177            pdv(:,:)    = pdv(:,:)    + zdvmr(:,:)
     1178            pdpsrf(:)   = pdpsrf(:)   + zdpsrfmr(:)
     1179            zdtsurf(:)  = zdtsurf(:)  + zdtsurfmr(:)
    11741180           
    11751181         endif
    … …  
    11791185  ! -----------------------------
    11801186
    1181         qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
     1187        qsurf(:,:) = qsurf(:,:) + ptimestep*dqsurf(:,:)
    11821188
    11831189         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
    … …  
    11951201      ! Increment surface temperature
    11961202
    1197       tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
     1203      tsurf(:)=tsurf(:)+ptimestep*zdtsurf(:)
    11981204
    11991205      ! Compute soil temperatures and subsurface heat flux.
    … …  
    12201226 
    12211227      ! Temperature, zonal and meridional winds.
    1222       zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
    1223       zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
    1224       zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
     1228      zt(:,:) = pt(:,:) + pdt(:,:)*ptimestep
     1229      zu(:,:) = pu(:,:) + pdu(:,:)*ptimestep
     1230      zv(:,:) = pv(:,:) + pdv(:,:)*ptimestep
    12251231
    12261232      ! Diagnostic.
    1227       zdtdyn(1:ngrid,1:nlayer)     = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
    1228       ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
    1229 
    1230       zdudyn(1:ngrid,1:nlayer)     = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep
    1231       zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer)
     1233      zdtdyn(:,:)     = (pt(:,:)-ztprevious(:,:)) / ptimestep
     1234      ztprevious(:,:) = zt(:,:)
     1235
     1236      zdudyn(:,:)     = (pu(:,:)-zuprevious(:,:)) / ptimestep
     1237      zuprevious(:,:) = zu(:,:)
    12321238
    12331239      if(firstcall)then
    1234          zdtdyn(1:ngrid,1:nlayer)=0.0
    1235          zdudyn(1:ngrid,1:nlayer)=0.0
     1240         zdtdyn(:,:)=0.0
     1241         zdudyn(:,:)=0.0
    12361242      endif
    12371243
    … …  
    12421248
    12431249      ! Tracers.
    1244       zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
     1250      zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
    12451251
    12461252      ! Surface pressure.
    1247       ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
     1253      ps(:) = pplev(:,1) + pdpsrf(:)*ptimestep
    12481254
    12491255
Note: See TracChangeset for help on using the changeset viewer.