Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp.F90 =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp.F90 (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp.F90 (revision 1415) @@ -23,4 +23,5 @@ ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +!!#include "cosp_defs.h" MODULE MOD_COSP USE MOD_COSP_TYPES @@ -128,12 +129,12 @@ ! and reff_zero == .false. Reff use in lidar and set to 0 for radar endif - if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed - print *, '---------- COSP ERROR ------------' - print *, '' - print *, 'use_reff==.true. but Reff is always zero' - print *, '' - print *, '----------------------------------' - stop - endif +! if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed +! print *, '---------- COSP ERROR ------------' +! print *, '' +! print *, 'use_reff==.true. but Reff is always zero' +! print *, '' +! print *, '----------------------------------' +! stop +! endif if ((.not. gbx%use_reff) .and. (reff_zero)) then ! No Reff in radar. Default in lidar gbx%Reff = DEFAULT_LIDAR_REFF @@ -322,4 +323,5 @@ integer :: Niter ! Number of calls to cosp_simulator integer :: i,j,k + integer :: I_HYDRO real,dimension(:,:),pointer :: column_frac_out ! Array with one column of frac_out integer,parameter :: scops_debug=0 ! set to non-zero value to print out inputs for debugging in SCOPS @@ -330,4 +332,5 @@ real,dimension(:,:),allocatable :: frac_ls,prec_ls,frac_cv,prec_cv ! Cloud/Precipitation fraction in each model level ! Levels are from SURFACE to TOA + real,dimension(:,:),allocatable :: rho ! (Npoints, Nlevels). Atmospheric dens type(cosp_sghydro) :: sghydro ! Subgrid info for hydrometeors en each iteration @@ -378,6 +381,6 @@ do k=1,Nlevels,1 do i=1,Ncolumns,1 - if (sgx%frac_out (j,i,Nlevels+1-k) .eq. 1) frac_ls(j,k)=frac_ls(j,k)+1. - if (sgx%frac_out (j,i,Nlevels+1-k) .eq. 2) frac_cv(j,k)=frac_cv(j,k)+1. + if (sgx%frac_out (j,i,Nlevels+1-k) == I_LSC) frac_ls(j,k)=frac_ls(j,k)+1. + if (sgx%frac_out (j,i,Nlevels+1-k) == I_CVC) frac_cv(j,k)=frac_cv(j,k)+1. if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 1) prec_ls(j,k)=prec_ls(j,k)+1. if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 2) prec_cv(j,k)=prec_cv(j,k)+1. @@ -414,5 +417,5 @@ !--------- Mixing ratios for clouds and Reff for Clouds and precip ------- column_frac_out => sgx%frac_out(:,k,:) - where (column_frac_out == 1) !+++++++++++ LS clouds ++++++++ + where (column_frac_out == I_LSC) !+++++++++++ LS clouds ++++++++ sghydro%mr_hydro(:,k,:,I_LSCLIQ) = gbx%mr_hydro(:,:,I_LSCLIQ) sghydro%mr_hydro(:,k,:,I_LSCICE) = gbx%mr_hydro(:,:,I_LSCICE) @@ -423,5 +426,5 @@ sghydro%Reff(:,k,:,I_LSSNOW) = gbx%Reff(:,:,I_LSSNOW) sghydro%Reff(:,k,:,I_LSGRPL) = gbx%Reff(:,:,I_LSGRPL) - elsewhere (column_frac_out == 2) !+++++++++++ CONV clouds ++++++++ + elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV clouds ++++++++ sghydro%mr_hydro(:,k,:,I_CVCLIQ) = gbx%mr_hydro(:,:,I_CVCLIQ) sghydro%mr_hydro(:,k,:,I_CVCICE) = gbx%mr_hydro(:,:,I_CVCICE) @@ -434,9 +437,9 @@ !--------- Precip ------- if (.not. gbx%use_precipitation_fluxes) then - where (column_frac_out == 1) !+++++++++++ LS Precipitation ++++++++ + where (column_frac_out == I_LSC) !+++++++++++ LS Precipitation ++++++++ sghydro%mr_hydro(:,k,:,I_LSRAIN) = gbx%mr_hydro(:,:,I_LSRAIN) sghydro%mr_hydro(:,k,:,I_LSSNOW) = gbx%mr_hydro(:,:,I_LSSNOW) sghydro%mr_hydro(:,k,:,I_LSGRPL) = gbx%mr_hydro(:,:,I_LSGRPL) - elsewhere (column_frac_out == 2) !+++++++++++ CONV Precipitation ++++++++ + elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV Precipitation ++++++++ sghydro%mr_hydro(:,k,:,I_CVRAIN) = gbx%mr_hydro(:,:,I_CVRAIN) sghydro%mr_hydro(:,k,:,I_CVSNOW) = gbx%mr_hydro(:,:,I_CVSNOW) @@ -488,5 +491,5 @@ sghydro%mr_hydro(:,:,:,I_LSRAIN),sghydro%mr_hydro(:,:,:,I_LSSNOW),sghydro%mr_hydro(:,:,:,I_LSGRPL), & sghydro%mr_hydro(:,:,:,I_CVRAIN),sghydro%mr_hydro(:,:,:,I_CVSNOW)) - endif + endif !++++++++++ CRM mode ++++++++++ else Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_constants.F90 =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_constants.F90 (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_constants.F90 (revision 1415) @@ -51,4 +51,10 @@ ! Number of possible output variables integer,parameter :: N_OUT_LIST = 27 + ! Value for forward model result from a level that is under the ground + real,parameter :: R_GROUND = -1.0E20 + + ! Stratiform and convective clouds in frac_out + integer, parameter :: I_LSC = 1, & ! Large-scale clouds + I_CVC = 2 ! Convective clouds !--- Radar constants @@ -56,5 +62,5 @@ integer,parameter :: DBZE_BINS = 15 ! Number of dBZe bins in histogram (cfad) real,parameter :: DBZE_MIN = -100.0 ! Minimum value for radar reflectivity - real,parameter :: DBZE_MAX = 30.0 ! Maximum value for radar reflectivity + real,parameter :: DBZE_MAX = 80.0 ! Maximum value for radar reflectivity real,parameter :: CFAD_ZE_MIN = -50.0 ! Lower value of the first CFAD Ze bin real,parameter :: CFAD_ZE_WIDTH = 5.0 ! Bin width (dBZe) @@ -69,6 +75,6 @@ integer,parameter :: LIDAR_NCAT = 4 integer,parameter :: PARASOL_NREFL = 5 ! parasol -! real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/0.0, 15.0, 30.0, 45.0, 60.0/) - real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/1.0, 2.0, 3.0, 4.0, 5.0/) + real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/0.0, 20.0, 40.0, 6.0, 80.0/) +! real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/1.0, 2.0, 3.0, 4.0, 5.0/) real,parameter :: DEFAULT_LIDAR_REFF = 30.0e-6 ! Default lidar effective radius Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_isccp_simulator.F90 =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_isccp_simulator.F90 (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_isccp_simulator.F90 (revision 1415) @@ -85,5 +85,5 @@ ! Change boxptop from hPa to Pa. This avoids using UDUNITS in CMOR - y%boxptop = y%boxptop*100.0 +! y%boxptop = y%boxptop*100.0 ! Check if there is any value slightly greater than 1 Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_simulator.F90 =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_simulator.F90 (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_simulator.F90 (revision 1415) @@ -65,20 +65,10 @@ !+++++++++ Radar model ++++++++++ if (cfg%Lradar_sim) then - call system_clock(t0,count_rate,count_max) call cosp_radar(gbx,sgx,sghydro,sgradar) - call system_clock(t1,count_rate,count_max) - print *, '%%%%%% Radar:', (t1-t0)*1.0/count_rate, ' s' - else - print *, '%%%%%% Radar not used' endif !+++++++++ Lidar model ++++++++++ if (cfg%Llidar_sim) then - call system_clock(t0,count_rate,count_max) call cosp_lidar(gbx,sgx,sghydro,sglidar) - call system_clock(t1,count_rate,count_max) - print *, '%%%%%% Lidar:', (t1-t0)*1.0/count_rate, ' s' - else - print *, '%%%%%% Lidar not used' endif @@ -86,41 +76,20 @@ !+++++++++ ISCCP simulator ++++++++++ if (cfg%Lisccp_sim) then - call system_clock(t0,count_rate,count_max) call cosp_isccp_simulator(gbx,sgx,isccp) - call system_clock(t1,count_rate,count_max) - print *, '%%%%%% ISCCP:', (t1-t0)*1.0/count_rate, ' s' - else - print *, '%%%%%% ISCCP not used' endif !+++++++++ MISR simulator ++++++++++ if (cfg%Lmisr_sim) then - call system_clock(t0,count_rate,count_max) call cosp_misr_simulator(gbx,sgx,misr) - call system_clock(t1,count_rate,count_max) - print *, '%%%%%% MISR:', (t1-t0)*1.0/count_rate, ' s' - else - print *, '%%%%%% MISR not used' endif !+++++++++++ Summary statistics +++++++++++ -! write(*,*) 'Stats:' -! read(*,*) c if (cfg%Lstats) then - call system_clock(t0,count_rate,count_max) call cosp_stats(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar) - call system_clock(t1,count_rate,count_max) - print *, '%%%%%% Stats:', (t1-t0)*1.0/count_rate, ' s' - endif - !+++++++++++ change of units after computation of statistics +++++++++++ - if (cfg%Llidar_sim) then - where((sglidar%beta_tot > 0.0) .and. (sglidar%beta_tot /= R_UNDEF)) - sglidar%beta_tot = log10(sglidar%beta_tot) - elsewhere - sglidar%beta_tot = R_UNDEF - end where +! print *, '%%%%%% Stats:', (t1-t0)*1.0/count_rate, ' s' endif + END SUBROUTINE COSP_SIMULATOR Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_stats.F90 =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_stats.F90 (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_stats.F90 (revision 1415) @@ -1,24 +1,24 @@ ! (c) British Crown Copyright 2008, the Met Office. ! All rights reserved. -! -! Redistribution and use in source and binary forms, with or without modification, are permitted +! +! Redistribution and use in source and binary forms, with or without modification, are permitted ! provided that the following conditions are met: -! -! * Redistributions of source code must retain the above copyright notice, this list +! +! * Redistributions of source code must retain the above copyright notice, this list ! of conditions and the following disclaimer. ! * Redistributions in binary form must reproduce the above copyright notice, this list -! of conditions and the following disclaimer in the documentation and/or other materials +! of conditions and the following disclaimer in the documentation and/or other materials ! provided with the distribution. -! * Neither the name of the Met Office nor the names of its contributors may be used -! to endorse or promote products derived from this software without specific prior written +! * Neither the name of the Met Office nor the names of its contributors may be used +! to endorse or promote products derived from this software without specific prior written ! permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR -! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR -! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER -! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @@ -29,7 +29,7 @@ ! Oct 2008 - J.-L. Dufresne - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS ! Oct 2008 - H. Chepfer - Added PARASOL reflectance arguments -! May 2010 - L. Fairhead - Optimisation of COSP_CHANGE_VERTICAL_GRID routine for NEC SX8 computer -! -! +! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations +! +! MODULE MOD_COSP_STATS USE MOD_COSP_CONSTANTS @@ -45,5 +45,5 @@ !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE COSP_STATS(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar) - + ! Input arguments type(cosp_gridbox),intent(in) :: gbx @@ -55,7 +55,7 @@ ! Output arguments type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar - type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar - - ! Local variables + type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar + + ! Local variables integer :: Npoints !# of grid points integer :: Nlevels !# of levels @@ -66,5 +66,5 @@ real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out real,dimension(:,:),allocatable :: ph_c,betamol_c - + Npoints = gbx%Npoints Nlevels = gbx%Nlevels @@ -73,5 +73,5 @@ Nlr = vgrid%Nlvgrid - if (cfg%Lcfad_lidarsr532) ok_lidar_cfad=.true. + if (cfg%Lcfad_Lidarsr532) ok_lidar_cfad=.true. if (vgrid%use_vgrid) then ! Statistics in a different vertical grid @@ -81,5 +81,4 @@ Ze_out = 0.0 betatot_out = 0.0 - betamol_in(:,1,:) = sglidar%beta_mol(:,:) betamol_out= 0.0 betamol_c = 0.0 @@ -89,5 +88,5 @@ !++++++++++++ Radar CFAD ++++++++++++++++ if (cfg%Lradar_sim) then - call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, & + call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, & Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.) stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, & @@ -96,4 +95,5 @@ !++++++++++++ Lidar CFAD ++++++++++++++++ if (cfg%Llidar_sim) then + betamol_in(:,1,:) = sglidar%beta_mol(:,:) call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,betamol_in, & Nlr,vgrid%zl,vgrid%zu,betamol_out) @@ -114,5 +114,5 @@ if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, & betatot_out,betamol_c,Ze_out, & - stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc) + stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc) ! Deallocate arrays at coarse resolution deallocate(Ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c) @@ -131,11 +131,11 @@ if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, & sglidar%beta_tot,sglidar%beta_mol,sgradar%Ze_tot, & - stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc) + stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc) endif ! Replace undef - where (stlidar%cfad_sr == LIDAR_UNDEF) stlidar%cfad_sr = R_UNDEF - where (stlidar%lidarcld == LIDAR_UNDEF) stlidar%lidarcld = R_UNDEF - where (stlidar%cldlayer == LIDAR_UNDEF) stlidar%cldlayer = R_UNDEF - where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF + where (stlidar%cfad_sr == LIDAR_UNDEF) stlidar%cfad_sr = R_UNDEF + where (stlidar%lidarcld == LIDAR_UNDEF) stlidar%lidarcld = R_UNDEF + where (stlidar%cldlayer == LIDAR_UNDEF) stlidar%cldlayer = R_UNDEF + where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF END SUBROUTINE COSP_STATS @@ -162,95 +162,93 @@ ! Local variables - integer :: i,j,k,c + integer :: i,j,k logical :: lunits - real :: ws(Npoints,M), ws_temp(Npoints,M) - real,dimension(Npoints,Nlevels) :: xl, xu ! Lower and upper boundaries of model grid - real,dimension(M) :: dz ! Layer depth - real,dimension(Npoints,Nlevels,M) :: w ! Weights to do the mean at each point - real,dimension(Npoints,Ncolumns,Nlevels) :: yp ! Variable to be changed to a different grid. - ! Local copy at a particular point. - ! This allows for change of units. - + + integer :: l + real,dimension(Npoints) :: ws,sumwyp + real,dimension(Npoints,Nlevels) :: xl,xu + real,dimension(Npoints,Nlevels) :: w + real,dimension(Npoints,Ncolumns,Nlevels) :: yp lunits=.false. if (present(log_units)) lunits=log_units - r(:,:,:) = 0.0 - yp(:,:,:) = y(:,:,:) - w(:,:,:) = 0.0 - ws(:,:) = 0.0 - ws_temp(:,:) = 0.0 - dz = zu - zl - + r(:,:,:) = R_GROUND + ! Vertical grid at that point + xl(:,:) = zhalf(:,:) + xu(:,1:Nlevels-1) = xl(:,2:Nlevels) + xu(:,Nlevels) = zfull(:,Nlevels) + zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric + yp(:,:,:) = y(:,:,:) ! Temporary variable to regrid ! Check for dBZ and change if necessary if (lunits) then where (y /= R_UNDEF) - yp = 10.0**(y/10.0) + yp = 10.0**(y/10.0) + elsewhere + yp = 0.0 end where endif - - ! Vertical grids - xl(:,:) = zhalf(:,:) - xu(:,1:Nlevels-1) = zhalf(:,2:Nlevels) - xu(:,Nlevels) = zfull(:,Nlevels) + zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric - ! Find weights - do k=1, M - do j=1, Nlevels - do i=1, Npoints + do k=1,M + ! Find weights + w(:,:) = 0.0 + do j=1,Nlevels + do i=1,Npoints if ((xl(i,j) < zl(k)).and.(xu(i,j) > zl(k)).and.(xu(i,j) <= zu(k))) then !xl(j)-----------------xu(j) ! zl(k)------------------------------zu(k) - w(i,j,k) = xu(i,j) - zl(k) + w(i,j) = xu(i,j) - zl(k) else if ((xl(i,j) >= zl(k)).and.(xu(i,j) <= zu(k))) then ! xl(j)-----------------xu(j) ! zl(k)------------------------------zu(k) - w(i,j,k) = xl(i,j+1) - xl(i,j) + w(i,j) = xu(i,j) - xl(i,j) else if ((xl(i,j) >= zl(k)).and.(xl(i,j) < zu(k)).and.(xu(i,j) >= zu(k))) then ! xl(j)-----------------xu(j) ! zl(k)------------------------------zu(k) - w(i,j,k) = zu(k) - xl(i,j) + w(i,j) = zu(k) - xl(i,j) else if ((xl(i,j) <= zl(k)).and.(xu(i,j) >= zu(k))) then ! xl(j)---------------------------xu(j) ! zl(k)--------------zu(k) - w(i,j,k) = dz(j) + w(i,j) = zu(k) - zl(k) + endif + enddo + enddo + ! Do the weighted mean + do j=1,Ncolumns + ws (:) = 0.0 + sumwyp(:) = 0.0 + do l=1,Nlevels + do i=1,Npoints + if (zu(k) > zhalf(i,1)) then ! Level above model bottom level + ws (i) = ws (i) + w(i,l) + sumwyp(i) = sumwyp(i) + w(i,l)*yp(i,j,l) + endif + enddo + enddo + do i=1,Npoints + if (zu(k) > zhalf(i,1)) then ! Level above model bottom level + if (ws(i) > 0.0) r(i,j,k) = sumwyp(i)/ws(i) endif enddo enddo enddo - - ! Do the weighted mean - do k=1,M - do i = 1, Npoints - ws(i,k) = sum(w(i,:,k)) - enddo - enddo - - ws_temp = 1. - where (ws(:,:) > 0.0) ws_temp(:,:)=ws(:,:) - - do c=1,Ncolumns - do k=1,M - do i = 1, Npoints - r(i,c,k) = sum(w(i,:,k)*yp(i,c,:))/ws_temp(i,k) - enddo - enddo - enddo - - do k=1,M - do i = 1, Npoints - if (ws(i,k) <= 0.0) r(i,:,k)=0.0 - enddo - enddo - ! Check for dBZ and change if necessary if (lunits) then - where (r(:,:,:) <= 0.0) - r(:,:,:) = R_UNDEF - elsewhere - r(:,:,:) = 10.0*log10(r(:,:,:)) - end where + do k=1,M + do j=1,Ncolumns + do i=1,Npoints + if (zu(k) > zhalf(i,1)) then ! Level above model bottom level + if (r(i,j,k) <= 0.0) then + r(i,j,k) = R_UNDEF + else + r(i,j,k) = 10.0*log10(r(i,j,k)) + endif + endif + enddo + enddo + enddo endif - -END SUBROUTINE COSP_CHANGE_VERTICAL_GRID + + + +END SUBROUTINE COSP_CHANGE_VERTICAL_GRID END MODULE MOD_COSP_STATS Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_utils.F90 =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_utils.F90 (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_utils.F90 (revision 1415) @@ -40,4 +40,46 @@ END INTERFACE CONTAINS + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +!------------------- SUBROUTINE COSP_PRECIP_MXRATIO -------------- +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +SUBROUTINE COSP_PRECIP_MXRATIO(Npoints,Nlevels,Ncolumns,p,T,prec_frac,prec_type, & + n_ax,n_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma1,gamma2, & + flux,mxratio) + + ! Input arguments, (IN) + integer,intent(in) :: Npoints,Nlevels,Ncolumns + real,intent(in),dimension(Npoints,Nlevels) :: p,T,flux + real,intent(in),dimension(Npoints,Ncolumns,Nlevels) :: prec_frac + real,intent(in) :: n_ax,n_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma1,gamma2,prec_type + ! Input arguments, (OUT) + real,intent(out),dimension(Npoints,Ncolumns,Nlevels) :: mxratio + ! Local variables + integer :: i,j,k + real :: sigma,one_over_xip1,xi,rho0,rho + + mxratio = 0.0 + + if (n_ax >= 0.0) then ! N_ax is used to control which hydrometeors need to be computed + !gamma1 = gamma(alpha_x + b_x + d_x + 1.0) + !gamma2 = gamma(alpha_x + b_x + 1.0) + xi = d_x/(alpha_x + b_x - n_bx + 1.0) + rho0 = 1.29 + sigma = (gamma2/(gamma1*c_x))*(n_ax*a_x*gamma2)**xi + one_over_xip1 = 1.0/(xi + 1.0) + + do k=1,Nlevels + do j=1,Ncolumns + do i=1,Npoints + if ((prec_frac(i,j,k)==prec_type).or.(prec_frac(i,j,k)==3.)) then + rho = p(i,k)/(287.05*T(i,k)) + mxratio(i,j,k)=(flux(i,k)*((rho/rho0)**g_x)*sigma)**one_over_xip1 + mxratio(i,j,k)=mxratio(i,j,k)/rho + endif + enddo + enddo + enddo + endif +END SUBROUTINE COSP_PRECIP_MXRATIO Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/icarus.F =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/icarus.F (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/icarus.F (revision 1415) @@ -33,5 +33,5 @@ &) -!Id: icarus.f,v 4.0 2009/02/12 13:59:20 hadmw Exp $ +!$Id: icarus.f,v 4.1 2010/05/27 16:30:18 hadmw Exp $ ! *****************************COPYRIGHT**************************** @@ -72,4 +72,5 @@ ! *****************************COPYRIGHT******************************* ! *****************************COPYRIGHT******************************* +! *****************************COPYRIGHT******************************* implicit none @@ -234,5 +235,4 @@ REAL attropmin (npoints) REAL atmax(npoints) - REAL atmin(npoints) REAL btcmin(npoints) REAL transmax(npoints) @@ -356,5 +356,4 @@ do j=1,npoints ptrop(j)=5000. - atmin(j) = 400. attropmin(j) = 400. atmax(j) = 0. @@ -373,8 +372,13 @@ itrop(j)=ilev end if - if (at(j,ilev) .gt. atmax(j)) atmax(j)=at(j,ilev) - if (at(j,ilev) .lt. atmin(j)) atmin(j)=at(j,ilev) enddo 12 continue + + do 13 ilev=1,nlev + do j=1,npoints + if (at(j,ilev) .gt. atmax(j) .and. + & ilev .ge. itrop(j)) atmax(j)=at(j,ilev) + enddo +13 continue end if Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/llnl_stats.F90 =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/llnl_stats.F90 (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/llnl_stats.F90 (revision 1415) @@ -24,4 +24,5 @@ MODULE MOD_LLNL_STATS + USE MOD_COSP_CONSTANTS IMPLICIT NONE @@ -62,6 +63,8 @@ do j = 1, Nlevels, 1 do k = 1, Ncolumns, 1 - do i = 1, Npoints, 1 - if ((x(i,k,j) >= xmin) .and. (x(i,k,j) <= xmax)) then + do i = 1, Npoints, 1 + if (x(i,k,j) == R_GROUND) then + cosp_cfad(i,:,j) = R_UNDEF + elseif ((x(i,k,j) >= xmin) .and. (x(i,k,j) <= xmax)) then ibin = ceiling((x(i,k,j) - bmin)/bwidth) if (ibin > Nbins) ibin = Nbins @@ -72,5 +75,5 @@ enddo !k enddo !j - cosp_cfad = cosp_cfad / Ncolumns + where ((cosp_cfad /= R_UNDEF).and.(cosp_cfad /= 0.0)) cosp_cfad = cosp_cfad / Ncolumns END FUNCTION COSP_CFAD @@ -98,6 +101,6 @@ integer :: pr,i,j -! lidar_only_freq_cloud = 0.0 -! tcc = 0.0 + lidar_only_freq_cloud = 0.0 + tcc = 0.0 do pr=1,Npoints do i=1,Ncolumns Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/lmd_ipsl_stats.F90 =================================================================== --- LMDZ4/branches/LMDZ4_AR5/libf/cosp/lmd_ipsl_stats.F90 (revision 1414) +++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/lmd_ipsl_stats.F90 (revision 1415) @@ -1,24 +1,24 @@ ! Copyright (c) 2009, Centre National de la Recherche Scientifique ! All rights reserved. -! -! Redistribution and use in source and binary forms, with or without modification, are permitted +! +! Redistribution and use in source and binary forms, with or without modification, are permitted ! provided that the following conditions are met: -! -! * Redistributions of source code must retain the above copyright notice, this list +! +! * Redistributions of source code must retain the above copyright notice, this list ! of conditions and the following disclaimer. ! * Redistributions in binary form must reproduce the above copyright notice, this list -! of conditions and the following disclaimer in the documentation and/or other materials +! of conditions and the following disclaimer in the documentation and/or other materials ! provided with the distribution. ! * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its -! contributors may be used to endorse or promote products derived from this software without +! contributors may be used to endorse or promote products derived from this software without ! specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR -! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND -! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR -! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER -! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @@ -39,13 +39,13 @@ ! ----------------------------------------------------------------------------------- ! Lidar outputs : -! +! ! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction ! from the lidar signals (ATB and molecular ATB) computed from model outputs ! + ! Compute CFADs of lidar scattering ratio SR and of depolarization index -! +! ! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France). ! -! December 2008, S. Bony, H. Chepfer and J-L. Dufresne : +! December 2008, S. Bony, H. Chepfer and J-L. Dufresne : ! - change of the cloud detection threshold S_cld from 3 to 5, for better ! with both day and night observations. The optical thinest clouds are missed. @@ -53,4 +53,10 @@ ! December 2008, A. Bodas-Salcedo: ! - Dimensions of pmol reduced to (npoints,llm) +! August 2009, A. Bodas-Salcedo: +! - Warning message regarding PARASOL being valid only over ocean deleted. +! February 2010, A. Bodas-Salcedo: +! - Undef passed into cosp_cfad_sr +! June 2010, T. Yokohata, T. Nishimura and K. Ogochi +! Optimisation of COSP_CFAD_SR ! ! Version 1.0 (June 2007) @@ -70,7 +76,7 @@ real undef ! undefined value - real pnorm(npoints,ncol,llm) ! lidar ATB + real pnorm(npoints,ncol,llm) ! lidar ATB real pmol(npoints,llm) ! molecular ATB - real land(npoints) ! Landmask [0 - Ocean, 1 - Land] + real land(npoints) ! Landmask [0 - Ocean, 1 - Land] real pplay(npoints,llm) ! pressure on model levels (Pa) logical ok_lidar_cfad ! true if lidar CFAD diagnostics need to be computed @@ -78,13 +84,13 @@ ! c outputs : - real lidarcld(npoints,llm) ! 3D "lidar" cloud fraction + real lidarcld(npoints,llm) ! 3D "lidar" cloud fraction real cldlayer(npoints,ncat) ! "lidar" cloud fraction (low, mid, high, total) - real cfad2(npoints,max_bin,llm) ! CFADs of SR - real srbval(max_bin) ! SR bins in CFADs + real cfad2(npoints,max_bin,llm) ! CFADs of SR + real srbval(max_bin) ! SR bins in CFADs real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance ! c threshold for cloud detection : - real S_clr - parameter (S_clr = 1.2) + real S_clr + parameter (S_clr = 1.2) real S_cld ! parameter (S_cld = 3.0) ! Previous thresold for cloud detection @@ -102,10 +108,5 @@ ! c 0- Initializations ! c ------------------------------------------------------- -! Parasol reflectance algorithm is not valid over land. Write -! a warning if there is no land. Landmask [0 - Ocean, 1 - Land] - IF ( MAXVAL(land(:)) .EQ. 0.0) THEN - WRITE (*,*) 'WARNING. PARASOL reflectance is not valid over land' & - & ,' and there is only land' - END IF +! ! Should be modified in future version @@ -116,5 +117,5 @@ ! c ------------------------------------------------------- ! -! where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) +! where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) ! x3d = pnorm/pmol ! elsewhere @@ -124,5 +125,5 @@ do ic = 1, ncol pnorm_c = pnorm(:,ic,:) - where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) + where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) x3d_c = pnorm_c/pmol elsewhere @@ -142,5 +143,5 @@ ! c ------------------------------------------------------- -! c 3- CFADs +! c 3- CFADs ! c ------------------------------------------------------- if (ok_lidar_cfad) then @@ -148,5 +149,5 @@ ! c CFADs of subgrid-scale lidar scattering ratios : ! c ------------------------------------------------------- - CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin, & + CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, & x3d, & S_att,S_clr,xmax,cfad2,srbval) @@ -172,16 +173,16 @@ ! if land=0 -> parasolrefl=parasolrefl parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) & - + (1.0 - MAX(1.0-land(:),0.0))*undef + + (1.0 - MAX(1.0-land(:),0.0))*undef enddo RETURN END SUBROUTINE diag_lidar - - + + !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !-------------------- FUNCTION COSP_CFAD_SR ------------------------ ! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins, & + SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, & x,S_att,S_clr,xmax,cfad,srbval) IMPLICIT NONE @@ -205,5 +206,5 @@ ! Input arguments integer Npoints,Ncolumns,Nlevels,Nbins - real xmax,S_att,S_clr + real xmax,S_att,S_clr,undef ! Input-output arguments real x(Npoints,Ncolumns,Nlevels) @@ -213,4 +214,5 @@ ! Local variables integer i, j, k, ib + real srbval_ext(0:Nbins) ! c ------------------------------------------------------- @@ -235,26 +237,26 @@ cfad(:,:,:) = 0.0 - + srbval_ext(1:Nbins) = srbval + srbval_ext(0) = -1.0 ! c ------------------------------------------------------- ! c c- Compute CFAD ! c ------------------------------------------------------- - do j = Nlevels, 1, -1 - do k = 1, Ncolumns - where ( x(:,k,j).le.srbval(1) ) & - cfad(:,1,j) = cfad(:,1,j) + 1.0 - enddo !k - enddo !j - - do ib = 2, Nbins - do j = Nlevels, 1, -1 - do k = 1, Ncolumns - where ( ( x(:,k,j).gt.srbval(ib-1) ) .and. ( x(:,k,j).le.srbval(ib) ) ) & - cfad(:,ib,j) = cfad(:,ib,j) + 1.0 - enddo !k - enddo !j - enddo !ib - - cfad(:,:,:) = cfad(:,:,:) / float(Ncolumns) + do j = 1, Nlevels + do ib = 1, Nbins + do k = 1, Ncolumns + do i = 1, Npoints + if (x(i,k,j) /= undef) then + if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) & + cfad(i,ib,j) = cfad(i,ib,j) + 1.0 + else + cfad(i,ib,j) = undef + endif + enddo + enddo + enddo + enddo + + where (cfad .ne. undef) cfad = cfad / float(Ncolumns) ! c ------------------------------------------------------- @@ -264,5 +266,5 @@ !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !-------------------- SUBROUTINE COSP_CLDFRAC ------------------- -! c Purpose: Cloud fraction diagnosed from lidar measurements +! c Purpose: Cloud fraction diagnosed from lidar measurements !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat, & @@ -288,7 +290,14 @@ real nsub(Npoints,Nlevels) - -! --------------------------------------------------------------- -! 1- initialization + real cldlay1(Npoints,Ncolumns) + real cldlay2(Npoints,Ncolumns) + real cldlay3(Npoints,Ncolumns) + real nsublay1(Npoints,Ncolumns) + real nsublay2(Npoints,Ncolumns) + real nsublay3(Npoints,Ncolumns) + + +! --------------------------------------------------------------- +! 1- initialization ! --------------------------------------------------------------- @@ -317,5 +326,5 @@ ! number of usefull sub-columns: - where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef) ) + where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef) ) srok(:,:,k)=1.0 elsewhere @@ -329,27 +338,45 @@ ! categories) : ! --------------------------------------------------------------- -! Test abderr + lidarcld = 0.0 + nsub = 0.0 + +!! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts. + cldlay1 = 0.0 + cldlay2 = 0.0 + cldlay3 = 0.0 + cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4 + nsublay1 = 0.0 + nsublay2 = 0.0 + nsublay3 = 0.0 + nsublay(:,:,4) = 0.0 do k = Nlevels, 1, -1 do ic = 1, Ncolumns do ip = 1, Npoints - iz=1 - p1 = pplay(ip,k) - if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds - iz=3 - else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid clouds - iz=2 + p1 = pplay(ip,k) + + if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds + cldlay3(ip,ic) = MAX(cldlay3(ip,ic), cldy(ip,ic,k)) + nsublay3(ip,ic) = MAX(nsublay3(ip,ic), srok(ip,ic,k)) + else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid clouds + cldlay2(ip,ic) = MAX(cldlay2(ip,ic), cldy(ip,ic,k)) + nsublay2(ip,ic) = MAX(nsublay2(ip,ic), srok(ip,ic,k)) + else + cldlay1(ip,ic) = MAX(cldlay1(ip,ic), cldy(ip,ic,k)) + nsublay1(ip,ic) = MAX(nsublay1(ip,ic), srok(ip,ic,k)) endif - cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k)) - cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k)) + cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k)) lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k) - - nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k)) nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k)) nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k) - enddo enddo enddo + cldlay(:,:,1) = cldlay1 + cldlay(:,:,2) = cldlay2 + cldlay(:,:,3) = cldlay3 + nsublay(:,:,1) = nsublay1 + nsublay(:,:,2) = nsublay2 + nsublay(:,:,3) = nsublay3 ! -- grid-box 3D cloud fraction @@ -369,6 +396,6 @@ do ic = 1, Ncolumns - cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz) - nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz) + cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz) + nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz) enddo @@ -383,4 +410,4 @@ END SUBROUTINE COSP_CLDFRAC ! --------------------------------------------------------------- - + END MODULE MOD_LMD_IPSL_STATS