Ignore:
Timestamp:
Aug 2, 2024, 2:12:03 PM (13 months ago)
Author:
abarral
Message:

Add missing klon on strataer_emiss_mod.F90
Correct various missing explicit declarations
Replace tabs by spaces (tabs are not part of the fortran charset)
Continue cleaning modules
Removed unused arguments and variables

Location:
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp
Files:
30 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/array_lib.F90

    r5099 r5158  
    141141  xsort = xarr(ist)
    142142 
    143   do i=1,nxx
     143  DO i=1,nxx
    144144    iloc = infind(xsort,xxarr(i),dist=d)
    145145    if (d > tol) then
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/congvec.h

    r5099 r5158  
    3737! *****************************COPYRIGHT*******************************
    3838
    39       do irand = 1, npoints
     39      DO irand = 1, npoints
    4040          ! Marsaglia CONG algorithm
    4141          seed(irand)=69069*seed(irand)+1234567
     
    4848      overflow_32=i2_16*i2_16
    4949      if ( overflow_32 .le. huge32 ) then
    50           do irand = 1, npoints
     50          DO irand = 1, npoints
    5151              ran(irand)=ran(irand)+1
    5252              ran(irand)=(ran(irand))-int(ran(irand))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/cosp_output_mod.F90

    r5157 r5158  
    7272         "clmcalipsoice", "CALIPSO Ice-Phase Mid Level Cloud Fraction", "%", (/ ('', i=1, 3) /))
    7373  TYPE(ctrl_outcosp), SAVE :: o_clmcalipsoliq = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    74          "clmcalipsoliq", "CALIPSO Liq-Phase Mid Level Cloud Fraction", "%", (/ ('', i=1, 3) /))                 
     74         "clmcalipsoliq", "CALIPSO Liq-Phase Mid Level Cloud Fraction", "%", (/ ('', i=1, 3) /))
    7575  TYPE(ctrl_outcosp), SAVE :: o_clhcalipsoice = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    7676         "clhcalipsoice", "CALIPSO Ice-Phase High Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 
     
    8080         "cltcalipsoice", "CALIPSO Ice-Phase Tot Level Cloud Fraction", "%", (/ ('', i=1, 3) /))
    8181  TYPE(ctrl_outcosp), SAVE :: o_cltcalipsoliq = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    82          "cltcalipsoliq", "CALIPSO Liq-Phase Tot Level Cloud Fraction", "%", (/ ('', i=1, 3) /))                 
     82         "cltcalipsoliq", "CALIPSO Liq-Phase Tot Level Cloud Fraction", "%", (/ ('', i=1, 3) /))
    8383  TYPE(ctrl_outcosp), SAVE :: o_cllcalipsoun = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    8484         "cllcalipsoun", "CALIPSO Undefined-Phase Low Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 
     
    9494         "clcalipsoliq", "Lidar Liq-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /))
    9595  TYPE(ctrl_outcosp), SAVE :: o_clcalipsoun = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    96          "clcalipsoun", "Lidar Undef-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /))     
     96         "clcalipsoun", "Lidar Undef-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /))
    9797  TYPE(ctrl_outcosp), SAVE :: o_clcalipsotmpice = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    9898         "clcalipsotmpice", "Lidar Ice-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /))
     
    115115         "clcalipsothin", "Lidar Thin profile Cloud Fraction", "%", (/ ('', i=1, 3) /))         !OPAQ
    116116  TYPE(ctrl_outcosp), SAVE :: o_clcalipsozopaque = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & !OPAQ
    117          "clcalipsozopaque", "Lidar z_opaque Fraction", "%", (/ ('', i=1, 3) /))                !OPAQ
     117         "clcalipsozopaque", "Lidar z_opaque Fraction", "%", (/ ('', i=1, 3) /))            !OPAQ
    118118  TYPE(ctrl_outcosp), SAVE :: o_clcalipsoopacity = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & !OPAQ
    119          "clcalipsoopacity", "Lidar opacity Fraction", "%", (/ ('', i=1, 3) /))                 !OPAQ
     119         "clcalipsoopacity", "Lidar opacity Fraction", "%", (/ ('', i=1, 3) /))                    !OPAQ
    120120
    121121  TYPE(ctrl_outcosp), SAVE :: o_proftemp = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &         !TIBO
     
    243243  integer                  :: Nlevlmdz, Ncolumns      ! Number of levels
    244244  real,dimension(Nlevlmdz) :: presnivs
    245   real                     :: dtime, freq_cosp, ecrit_day, ecrit_hf, ecrit_mth
     245  REAL                     :: dtime, freq_cosp, ecrit_day, ecrit_hf, ecrit_mth
    246246  logical                  :: ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, use_vgrid, ok_all_xml                   
    247247  type(cosp_vgrid)   :: vgrid   ! Information on vertical grid of stats
     
    250250!!! Variables locales
    251251  integer                   :: idayref, iff, ii
    252   real                      :: zjulian,zjulian_start
     252  REAL                      :: zjulian,zjulian_start
    253253  real,dimension(Ncolumns)  :: column_ax
    254254  real,dimension(DBZE_BINS) ::  dbze_ax
     
    271271
    272272!! Definition valeurs axes
    273     do ii=1,Ncolumns
     273    DO ii=1,Ncolumns
    274274      column_ax(ii) = real(ii)
    275275    enddo
    276276
    277     do i=1,DBZE_BINS
     277    DO i=1,DBZE_BINS
    278278     dbze_ax(i) = CFAD_ZE_MIN + CFAD_ZE_WIDTH*(i - 0.5)
    279279    enddo
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/cosp_output_write_mod.F90

    r5157 r5158  
    3131!!! Variables d'entree
    3232  integer               :: itap, Nlevlmdz, Ncolumns, Npoints
    33   real                  :: freq_COSP, dtime, missing_val, missing_cosp
     33  REAL                  :: freq_COSP, dtime, missing_val, missing_cosp
    3434  type(cosp_config)     :: cfg     ! Control outputs
    3535  type(cosp_gridbox)    :: gbx     ! Gridbox information. Input for COSP
     
    161161
    162162   IF (using_xios) THEN
    163      do icl=1,SR_BINS
     163     DO icl=1,SR_BINS
    164164        tmp_fi4da_cfadL(:,:,icl)=stlidar%cfad_sr(:,icl,:)
    165165     enddo
     
    169169   ELSE
    170170     if (cfg%LcfadLidarsr532) then
    171        do icl=1,SR_BINS
     171       DO icl=1,SR_BINS
    172172          CALL histwrite3d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr(:,icl,:),nvert,icl)
    173173       enddo
    174174     endif
    175175     if (cfg%LprofSR) then
    176        do icl=1,Ncolumns                                                              !TIBO
     176       DO icl=1,Ncolumns                                                              !TIBO
    177177          CALL histwrite3d_cosp(o_profSR,stlidar%profSR(:,icl,:),nvert,icl)           !TIBO
    178178       enddo                                                                          !TIBO
     
    183183
    184184  if (cfg%LparasolRefl) then
    185     do k=1,PARASOL_NREFL
    186      do ip=1, Npoints
     185    DO k=1,PARASOL_NREFL
     186     DO ip=1, Npoints
    187187      if (stlidar%cldlayer(ip,4).gt.0.01.and.stlidar%parasolrefl(ip,k).ne.missing_val) then
    188188        parasolcrefl(ip,k)=(stlidar%parasolrefl(ip,k)-0.03*(1.-stlidar%cldlayer(ip,4)))/ &
     
    203203   ELSE
    204204     if (cfg%Latb532) then
    205        do icl=1,Ncolumns
     205       DO icl=1,Ncolumns
    206206          CALL histwrite3d_cosp(o_atb532,sglidar%beta_tot(:,icl,:),nvertmcosp,icl)
    207207       enddo
     
    218218   where(stradar%cfad_ze == R_UNDEF) stradar%cfad_ze = missing_val
    219219   IF (using_xios) THEN
    220      do icl=1,DBZE_BINS
     220     DO icl=1,DBZE_BINS
    221221       tmp_fi4da_cfadR(:,:,icl)=stradar%cfad_ze(:,icl,:)
    222222     enddo
     
    226226   ELSE
    227227     if (cfg%Ldbze94) then
    228        do icl=1,Ncolumns
     228       DO icl=1,Ncolumns
    229229         CALL histwrite3d_cosp(o_dbze94,sgradar%Ze_tot(:,icl,:),nvert,icl)
    230230       enddo
    231231     endif
    232232     if (cfg%LcfadDbze94) then
    233        do icl=1,DBZE_BINS
     233       DO icl=1,DBZE_BINS
    234234         CALL histwrite3d_cosp(o_cfadDbze94,stradar%cfad_ze(:,icl,:),nvert,icl)
    235235       enddo
     
    266266   ELSE
    267267     if (cfg%Lclisccp) then
    268        do icl=1,7
     268       DO icl=1,7
    269269         CALL histwrite3d_cosp(o_clisccp2,isccp%fq_isccp(:,icl,:),nvertisccp,icl)
    270270       enddo
     
    287287
    288288   IF (using_xios) THEN
    289      do icl=1,MISR_N_CTH
     289     DO icl=1,MISR_N_CTH
    290290        tmp_fi4da_misr(:,icl,:)=misr%fq_MISR(:,:,icl)
    291291     enddo
     
    294294   ELSE
    295295     if (cfg%LclMISR) then
    296       do icl=1,7
     296      DO icl=1,7
    297297        CALL histwrite3d_cosp(o_clMISR,misr%fq_MISR(:,icl,:),nvertmisr,icl)
    298298      enddo
     
    367367   ELSE
    368368     if (cfg%Lclmodis) then
    369        do icl=1,7
     369       DO icl=1,7
    370370         CALL histwrite3d_cosp(o_clmodis, &
    371371         modis%Optical_Thickness_vs_Cloud_Top_Pressure(:,icl,:),nvertisccp,icl)
     
    385385   ELSE
    386386     if (cfg%Lcrimodis) then
    387        do icl=1,7
     387       DO icl=1,7
    388388         CALL histwrite3d_cosp(o_crimodis, &
    389389            modis%Optical_Thickness_vs_ReffIce(:,icl,:),nvertReffIce,icl)
     
    391391     endif
    392392     if (cfg%Lcrlmodis) then
    393        do icl=1,7
     393       DO icl=1,7
    394394         CALL histwrite3d_cosp(o_crlmodis, &
    395395            modis%Optical_Thickness_vs_ReffLiq(:,icl,:),nvertReffLiq,icl)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/cosp_read_otputkeys.F90

    r5133 r5158  
    2828
    2929               
    30    do i=1,N_OUT_LIST
     30   DO i=1,N_OUT_LIST
    3131      cfg%out_list(i)=''
    3232   enddo
     
    130130
    131131
    132    do i=1,N_OUT_LIST
     132   DO i=1,N_OUT_LIST
    133133      cfg%out_list(i)=''
    134134   enddo
     
    269269             LprofSR,Lproftemp                                                                        !TIBO (2)
    270270   
    271   do i=1,N_OUT_LIST
     271  DO i=1,N_OUT_LIST
    272272    cfg%out_list(i)=''
    273273  enddo
     
    771771  IF (using_xios) THEN
    772772   
    773     do i=1,N_OUT_LIST
     773    DO i=1,N_OUT_LIST
    774774      cfg%out_list(i)=''
    775775    enddo
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/dsd.F90

    r5099 r5158  
    311311      lidx = infind(D,dmin)
    312312      uidx = infind(D,dmax)   
    313       do k=lidx,uidx
     313      DO k=lidx,uidx
    314314 
    315315        N(k) = ( &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/format_input.F90

    r5099 r5158  
    5252
    5353  nprof = size(hgt_matrix,1)
    54   do i=1,nprof
     54  DO i=1,nprof
    5555    call lin_interpolate(env_t_matrix(i,:),env_hgt_matrix(i,:), &
    5656      t_matrix(i,:),hgt_matrix(i,:),1000._KR8)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/gases.F90

    r5099 r5158  
    126126  sumo = 0.
    127127  aux1 = 1.1*e_th
    128   do i=1,nbands_o2
     128  DO i=1,nbands_o2
    129129    aux2 = f/v0(i)
    130130    aux3 = v0(i)-f
     
    151151  sumo = 0.
    152152  aux1 = 4.8*e_th
    153   do i=1,nbands_h2o
     153  DO i=1,nbands_h2o
    154154    aux2 = f/v1(i)
    155155    aux3 = v1(i)-f
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/lidar_simulator.F90

    r5099 r5158  
    324324! altitude at half pressure levels:
    325325      zheight(:,1) = 0.0
    326       do k = 2, nlev+1
     326      DO k = 2, nlev+1
    327327        zheight(:,k) = zheight(:,k-1) &
    328328                  -(presf(:,k)-presf(:,k-1))/(rhoair(:,k-1)*9.81)
     
    341341
    342342! polynomes kp_lidar derived from Mie theory:
    343       do i = 1, npart
     343      DO i = 1, npart
    344344       where ( rad_part(:,:,i).gt.0.0)
    345345         kp_part(:,:,i) = &
     
    361361
    362362! alpha of particles in each subcolumn:
    363       do i = 1, npart
     363      DO i = 1, npart
    364364        where ( rad_part(:,:,i).gt.0.0)
    365365          alpha_part(:,:,i) = 3.0/4.0 * Qscat &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/math_lib.F90

    r5099 r5158  
    4949      ga=1.0
    5050      m1=x-1
    51       do k=2,m1
     51      DO k=2,m1
    5252        ga=ga*k
    5353      enddo
     
    6060      m=int(z)
    6161      r=1.0
    62       do k=1,m
     62      DO k=1,m
    6363        r=r*(z-k)
    6464      enddo
     
    8080          -.206d-13, -.54d-14, .14d-14, .1d-15/
    8181    gr=g(26)
    82     do k=25,1,-1
     82    DO k=25,1,-1
    8383      gr=gr*z+g(k)
    8484    enddo
     
    146146  else
    147147     sumo = 0.
    148      do j=i1,i2
     148     DO j=i1,i2
    149149       deltah = abs(s(i1+1)-s(i1))
    150150       sumo = sumo + f(j)*deltah
     
    268268  end if
    269269 
    270   do i = 2, ntab
     270  DO i = 2, ntab
    271271    if ( xtab(i) <= xtab(i-1) ) then
    272272       lerror = .true.
     
    311311  ihi = ntab
    312312
    313   do i = 1, ntab
     313  DO i = 1, ntab
    314314    if ( a <= xtab(i) ) then
    315315      exit
     
    321321  ilo = min ( ilo, ntab - 1 )
    322322
    323   do i = 1, ntab
     323  DO i = 1, ntab
    324324    if ( xtab(i) <= b ) then
    325325      exit
     
    335335  sum1 = 0.0D+00
    336336 
    337   do i = ilo, ihi
     337  DO i = ilo, ihi
    338338 
    339339    x1 = xtab(i-1)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp.F90

    r5099 r5158  
    6868  integer,dimension(2) :: ix,iy
    6969  logical :: reff_zero
    70   real :: maxp,minp
     70  REAL :: maxp,minp
    7171  integer,dimension(:),allocatable :: & ! Dimensions nPoints
    7272                  seed    !  It is recommended that the seed is set to a different value for each model
     
    9595
    9696!++++++++++ Depth of model layers ++++++++++++
    97   do i=1,Nlevels-1
     97  DO i=1,Nlevels-1
    9898    gbx%dlev(:,i) = gbx%zlev_half(:,i+1) - gbx%zlev_half(:,i)
    9999  enddo
     
    195195        Niter = gbx%Npoints/gbx%Npoints_it ! Integer division
    196196        if (Niter*gbx%Npoints_it < gbx%Npoints) Niter = Niter + 1
    197         do i=1,Niter
     197        DO i=1,Niter
    198198            i_first = (i-1)*gbx%Npoints_it + 1
    199199            i_last  = i_first + gbx%Npoints_it - 1
     
    409409       
    410410        ! Precipitation fraction
    411         do j=1,Npoints,1
    412         do k=1,Nlevels,1
    413             do i=1,Ncolumns,1
     411        DO j=1,Npoints,1
     412        DO k=1,Nlevels,1
     413            DO i=1,Ncolumns,1
    414414                if (sgx%frac_out (j,i,Nlevels+1-k) == I_LSC) frac_ls(j,k)=frac_ls(j,k)+1.
    415415                if (sgx%frac_out (j,i,Nlevels+1-k) == I_CVC) frac_cv(j,k)=frac_cv(j,k)+1.
     
    434434        else
    435435            ! This is done within a loop (unvectorized) over nPoints to save memory
    436             do j=1,Npoints
     436            DO j=1,Npoints
    437437                sgx%frac_out(j,:,1:Nlevels)  = sgx%frac_out(j,:,Nlevels:1:-1)
    438438                sgx%prec_frac(j,:,1:Nlevels) = sgx%prec_frac(j,:,Nlevels:1:-1)
     
    445445        ! Populate the subgrid arrays
    446446        call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,Nhydro,sghydro)
    447         do k=1,Ncolumns
     447        DO k=1,Ncolumns
    448448            !--------- Mixing ratios for clouds and Reff for Clouds and precip -------
    449449            column_frac_out => sgx%frac_out(:,k,:)
     
    498498        enddo
    499499        ! convert the mixing ratio and precipitation flux from gridbox mean to the fraction-based values
    500         do k=1,Nlevels
    501             do j=1,Npoints
     500        DO k=1,Nlevels
     501            DO j=1,Npoints
    502502                !--------- Clouds -------
    503503                if (frac_ls(j,k) .ne. 0.) then
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_constants.F90

    r5099 r5158  
    164164    integer :: HCLASS_TYPE(N_HYDRO),HCLASS_PHASE(N_HYDRO)
    165165
    166     real :: HCLASS_DMIN(N_HYDRO),HCLASS_DMAX(N_HYDRO), &
     166    REAL :: HCLASS_DMIN(N_HYDRO),HCLASS_DMAX(N_HYDRO), &
    167167            HCLASS_APM(N_HYDRO),HCLASS_BPM(N_HYDRO),HCLASS_RHO(N_HYDRO), &
    168168            HCLASS_P1(N_HYDRO),HCLASS_P2(N_HYDRO),HCLASS_P3(N_HYDRO)
     
    274274    integer :: HCLASS_TYPE(N_HYDRO),HCLASS_PHASE(N_HYDRO)
    275275
    276     real :: HCLASS_DMIN(N_HYDRO),HCLASS_DMAX(N_HYDRO), &           
     276    REAL :: HCLASS_DMIN(N_HYDRO),HCLASS_DMAX(N_HYDRO), &
    277277            HCLASS_APM(N_HYDRO),HCLASS_BPM(N_HYDRO),HCLASS_RHO(N_HYDRO), &
    278278            HCLASS_P1(N_HYDRO),HCLASS_P2(N_HYDRO),HCLASS_P3(N_HYDRO)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_isccp_simulator.F90

    r5099 r5158  
    4545  ! Local variables
    4646  integer :: Nlevels,Npoints
    47   real :: pfull(gbx%Npoints, gbx%Nlevels)
    48   real :: phalf(gbx%Npoints, gbx%Nlevels + 1)
    49   real :: qv(gbx%Npoints, gbx%Nlevels)
    50   real :: cc(gbx%Npoints, gbx%Nlevels)
    51   real :: conv(gbx%Npoints, gbx%Nlevels)
    52   real :: dtau_s(gbx%Npoints, gbx%Nlevels)
    53   real :: dtau_c(gbx%Npoints, gbx%Nlevels)
    54   real :: at(gbx%Npoints, gbx%Nlevels)
    55   real :: dem_s(gbx%Npoints, gbx%Nlevels)
    56   real :: dem_c(gbx%Npoints, gbx%Nlevels)
    57   real :: frac_out(gbx%Npoints, gbx%Ncolumns, gbx%Nlevels)
     47  REAL :: pfull(gbx%Npoints, gbx%Nlevels)
     48  REAL :: phalf(gbx%Npoints, gbx%Nlevels + 1)
     49  REAL :: qv(gbx%Npoints, gbx%Nlevels)
     50  REAL :: cc(gbx%Npoints, gbx%Nlevels)
     51  REAL :: conv(gbx%Npoints, gbx%Nlevels)
     52  REAL :: dtau_s(gbx%Npoints, gbx%Nlevels)
     53  REAL :: dtau_c(gbx%Npoints, gbx%Nlevels)
     54  REAL :: at(gbx%Npoints, gbx%Nlevels)
     55  REAL :: dem_s(gbx%Npoints, gbx%Nlevels)
     56  REAL :: dem_c(gbx%Npoints, gbx%Nlevels)
     57  REAL :: frac_out(gbx%Npoints, gbx%Ncolumns, gbx%Nlevels)
    5858  integer :: sunlit(gbx%Npoints)
    5959 
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_lidar.F90

    r5099 r5158  
    5454  ! Local variables
    5555  integer :: i
    56   real :: presf(sgx%Npoints, sgx%Nlevels + 1)
     56  REAL :: presf(sgx%Npoints, sgx%Nlevels + 1)
    5757  real,dimension(sgx%Npoints, sgx%Nlevels) :: lsca,mr_ll,mr_li,mr_cl,mr_ci
    5858  real,dimension(sgx%Npoints, sgx%Nlevels) :: beta_tot,tau_tot
     
    6363  presf(:,sgx%Nlevels + 1) = 0.0
    6464  lsca = gbx%tca-gbx%cca
    65   do i=1,sgx%Ncolumns
     65  DO i=1,sgx%Ncolumns
    6666      ! Temporary arrays for simulator call
    6767      mr_ll(:,:) = sghydro%mr_hydro(:,i,:,I_LSCLIQ)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_misr_simulator.F90

    r5099 r5158  
    4949  ! Local variables
    5050  integer :: Nlevels,Npoints
    51   real :: dtau_s(gbx%Npoints, gbx%Nlevels)
    52   real :: dtau_c(gbx%Npoints, gbx%Nlevels)
    53   real :: at(gbx%Npoints, gbx%Nlevels)
    54   real :: frac_out(gbx%Npoints, gbx%Ncolumns, gbx%Nlevels)
     51  REAL :: dtau_s(gbx%Npoints, gbx%Nlevels)
     52  REAL :: dtau_c(gbx%Npoints, gbx%Nlevels)
     53  REAL :: at(gbx%Npoints, gbx%Nlevels)
     54  REAL :: frac_out(gbx%Npoints, gbx%Ncolumns, gbx%Nlevels)
    5555  integer :: sunlit(gbx%Npoints)
    5656 
    57   real :: zfull(gbx%Npoints, gbx%Nlevels) !  height (in meters) of full model levels (i.e. midpoints)
     57  REAL :: zfull(gbx%Npoints, gbx%Nlevels) !  height (in meters) of full model levels (i.e. midpoints)
    5858                                          !  zfull(npoints,1)    is    top level of model
    5959                                          !  zfull(npoints,nlev) is bottom level of model
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_modis_simulator.F90

    r5099 r5158  
    162162
    163163      ! Loop version of spread above - intrinsic doesn't work on certain platforms.
    164       do k = 1, nLevels
    165         do j = 1, nSubCols
    166           do i = 1, nSunlit
     164      DO k = 1, nLevels
     165        DO j = 1, nSubCols
     166          DO i = 1, nSunlit
    167167            if(subCols%frac_out(sunlit(i), j, k) == I_LSC) then
    168168              opticalThickness(i, j, k) = gridBox%dtau_s(sunlit(i), k)
     
    186186
    187187      ! Loop version of spread above - intrinsic doesn't work on certain platforms.
    188       do k = 1, nLevels
    189         do j = 1, nSubCols
    190           do i = 1, nSunlit
     188      DO k = 1, nLevels
     189        DO j = 1, nSubCols
     190          DO i = 1, nSunlit
    191191            if(subCols%frac_out(sunlit(i), j, k) == I_CVC) opticalThickness(i, j, k) = gridBox%dtau_c(sunlit(i), k)
    192192          end do
     
    205205      isccpCloudTopPressure(:, :) = isccpSim%boxptop(sunlit(:), :)
    206206     
    207       do i = 1, nSunlit
     207      DO i = 1, nSunlit
    208208        call modis_L2_simulator(temperature(i, :), pressureLayers(i, :), pressureLevels(i, :),     &
    209209                                opticalThickness(i, :, :), cloudWater(i, :, :), cloudIce(i, :, :), &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_radar.F90

    r5099 r5158  
    139139
    140140  ! set flag denoting position of radar relative to hgt_matrix orientation
    141           ngate = size(hgt_matrix,2)
     141      ngate = size(hgt_matrix,2)
    142142
    143           hgt_descending = hgt_matrix(1,1) > hgt_matrix(1,ngate)
     143      hgt_descending = hgt_matrix(1,1) > hgt_matrix(1,ngate)
    144144
    145           if ( &
    146              (gbx%surface_radar == 1 .and. hgt_descending) .or.  &
    147              (gbx%surface_radar == 0 .and. (.not. hgt_descending)) &
    148              ) &
    149           then
    150             gbx%hp%radar_at_layer_one = .false.
    151           else
    152             gbx%hp%radar_at_layer_one = .true.
    153           endif
     145      if ( &
     146         (gbx%surface_radar == 1 .and. hgt_descending) .or.  &
     147         (gbx%surface_radar == 0 .and. (.not. hgt_descending)) &
     148         ) &
     149      then
     150        gbx%hp%radar_at_layer_one = .false.
     151      else
     152        gbx%hp%radar_at_layer_one = .true.
     153      endif
    154154
    155155  ! ----- loop over subcolumns -----
    156   do pr=1,sgx%Ncolumns
     156  DO pr=1,sgx%Ncolumns
    157157
    158158      !  NOTE:
     
    160160      !  only hydrometeor profiles will be different for each subgridbox
    161161
    162          do i=1,gbx%Nhydro
     162         DO i=1,gbx%Nhydro
    163163            hm_matrix(i,:,:) = sghydro%mr_hydro(:,pr,:,i)*1000.0 ! Units from kg/kg to g/kg
    164164            if (gbx%use_reff) then
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_stats.F90

    r5099 r5158  
    203203   logical :: lunits
    204204   integer :: l
    205    real :: w ! Weight
    206    real :: dbb, dtb, dbt, dtt ! Distances between edges of both grids
     205   REAL :: w ! Weight
     206   REAL :: dbb, dtb, dbt, dtt ! Distances between edges of both grids
    207207   integer :: Nw  ! Number of weights
    208    real :: wt  ! Sum of weights
     208   REAL :: wt  ! Sum of weights
    209209   real,dimension(Nlevels) :: oldgrid_bot,oldgrid_top ! Lower and upper boundaries of model grid
    210    real :: yp ! Local copy of y at a particular point.
     210   REAL :: yp ! Local copy of y at a particular point.
    211211              ! This allows for change of units.
    212212
     
    216216   r = 0.0
    217217
    218    do i=1,Npoints
     218   DO i=1,Npoints
    219219     ! Calculate tops and bottoms of new and old grids
    220220     oldgrid_bot = zhalf(i,:)
     
    223223     l = 0 ! Index of level in the old grid
    224224     ! Loop over levels in the new grid
    225      do k = 1,Nglevels
     225     DO k = 1,Nglevels
    226226       Nw = 0 ! Number of weigths
    227227       wt = 0.0 ! Sum of weights
    228228       ! Loop over levels in the old grid and accumulate total for weighted average
    229        do
     229       DO
    230230         l = l + 1
    231231         w = 0.0 ! Initialise weight to 0
     
    254254             Nw = Nw + 1
    255255             wt = wt + w
    256              do j=1,Ncolumns
     256             DO j=1,Ncolumns
    257257               if (lunits) then
    258258                 if (y(i,j,l) /= R_UNDEF) then
     
    273273       ! Calculate average in new grid
    274274       if (Nw > 0) then
    275          do j=1,Ncolumns
     275         DO j=1,Ncolumns
    276276           r(i,j,k) = r(i,j,k)/wt
    277277         enddo
     
    281281
    282282   ! Set points under surface to R_UNDEF, and change to dBZ if necessary
    283    do k=1,Nglevels
    284      do j=1,Ncolumns
    285        do i=1,Npoints
     283   DO k=1,Nglevels
     284     DO j=1,Ncolumns
     285       DO i=1,Npoints
    286286         if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level
    287287           if (lunits) then
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_types.F90

    r5099 r5158  
    4444                Lclcalipsoliq,Lclcalipsoice,Lclcalipsoun, &
    4545                Lclcalipsotmp,Lclcalipsotmpliq,Lclcalipsotmpice,Lclcalipsotmpun, &
    46                       Lcltcalipsoliq,Lcltcalipsoice,Lcltcalipsoun, &
     46                  Lcltcalipsoliq,Lcltcalipsoice,Lcltcalipsoun, &
    4747                Lclhcalipsoliq,Lclhcalipsoice,Lclhcalipsoun, &
    4848                Lclmcalipsoliq,Lclmcalipsoice,Lclmcalipsoun, &
     
    274274   
    275275    ! Radar ancillary info
    276     real :: radar_freq, & ! Radar frequency [GHz]
     276    REAL :: radar_freq, & ! Radar frequency [GHz]
    277277            k2 ! |K|^2, -1=use frequency dependent default
    278278    integer :: surface_radar, & ! surface=1, spaceborne=0
     
    383383                                 ! 2 = experimental setting 
    384384    integer :: isccp_overlap !  overlap type (1=max, 2=rand, 3=max/rand)
    385     real :: isccp_emsfc_lw      ! 10.5 micron emissivity of surface (fraction)
     385    REAL :: isccp_emsfc_lw      ! 10.5 micron emissivity of surface (fraction)
    386386 
    387387    ! RTTOV inputs/options
     
    392392    integer, dimension(:), pointer :: Ichan   ! Channel numbers
    393393    real,    dimension(:), pointer :: Surfem  ! Surface emissivity
    394     real    :: ZenAng ! Satellite Zenith Angles
    395     real :: co2,ch4,n2o,co ! Mixing ratios of trace gases
     394    REAL    :: ZenAng ! Satellite Zenith Angles
     395    REAL :: co2,ch4,n2o,co ! Mixing ratios of trace gases
    396396
    397397  END TYPE COSP_GRIDBOX
     
    550550    ! Local variables
    551551    integer :: i
    552     real :: zstep
     552    REAL :: zstep
    553553   
    554554    x%use_vgrid  = use_vgrid
     
    587587         zstep = 20000.0/x%Nlvgrid
    588588      endif
    589       do i=1,x%Nlvgrid
     589      DO i=1,x%Nlvgrid
    590590         x%zl(i) = (i-1)*zstep
    591591         x%zu(i) = i*zstep
     
    654654    x%tau_tot    = 0.0
    655655    x%refl       = 0.0 ! parasol
    656     x%temp_tot          = 0.0
    657     x%betaperp_tot      = 0.0   
     656    x%temp_tot       = 0.0
     657    x%betaperp_tot     = 0.0
    658658  END SUBROUTINE CONSTRUCT_COSP_SGLIDAR
    659659
     
    12041204    ! y%hp%idd     = x%hp%idd
    12051205    sz = shape(x%hp%Z_scale_flag)
    1206     do k=1,sz(3)
    1207       do j=1,sz(2)
    1208         do i=1,sz(1)
     1206    DO k=1,sz(3)
     1207      DO j=1,sz(2)
     1208        DO i=1,sz(1)
    12091209           if (x%hp%N_scale_flag(i,k))   y%hp%N_scale_flag(i,k)      = .true.
    12101210           if (x%hp%Z_scale_flag(i,j,k)) y%hp%Z_scale_flag(i,j,k)    = .true.
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_utils.F90

    r5099 r5158  
    5858    ! Local variables
    5959    integer :: i,j,k
    60     real :: sigma,one_over_xip1,xi,rho0,rho,lambda_x,gamma_4_3_2,delta
    61     real :: seuil
     60    REAL :: sigma,one_over_xip1,xi,rho0,rho,lambda_x,gamma_4_3_2,delta
     61    REAL :: seuil
    6262
    6363    if (ok_debug_cosp) then
     
    7777        delta = (alpha_x + b_x + d_x - n_bx + 1.0)
    7878       
    79         do k=1,Nlevels
    80             do j=1,Ncolumns
    81                 do i=1,Npoints
     79        DO k=1,Nlevels
     80            DO j=1,Ncolumns
     81                DO i=1,Npoints
    8282                    if ((prec_frac(i,j,k)==prec_type).or.(prec_frac(i,j,k)==3.)) then
    8383                        rho = p(i,k)/(287.05*T(i,k))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_llnl_stats.F90

    r5099 r5158  
    6767   ! Valid x values smaller than bmin and larger than bmax are set
    6868   ! into the smallest bin and largest bin, respectively.
    69    do j = 1, Nlevels, 1
    70       do k = 1, Ncolumns, 1
    71          do i = 1, Npoints, 1
     69   DO j = 1, Nlevels, 1
     70      DO k = 1, Ncolumns, 1
     71         DO i = 1, Npoints, 1
    7272            if (x(i,k,j) == R_GROUND) then
    7373               cosp_cfad(i,:,j) = R_UNDEF
     
    101101
    102102   ! local variables
    103    real :: sc_ratio
    104    real :: s_cld, s_att
     103   REAL :: sc_ratio
     104   REAL :: s_cld, s_att
    105105   parameter (S_cld = 5.0)
    106106   parameter (s_att = 0.01)
     
    111111   lidar_only_freq_cloud = 0.0
    112112   tcc = 0.0
    113    do pr=1,Npoints
    114      do i=1,Ncolumns
     113   DO pr=1,Npoints
     114     DO i=1,Ncolumns
    115115       flag_sat = 0
    116116       flag_cld = 0
    117        do j=Nlevels,1,-1 !top->surf
     117       DO j=Nlevels,1,-1 !top->surf
    118118        sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j)
    119119        if ((sc_ratio .le. s_att) .and. (flag_sat .eq. 0)) flag_sat = j
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_lmd_ipsl_stats.F90

    r5099 r5158  
    131131
    132132! SR detection threshold
    133       real, parameter  ::  S_cld_att = 30. ! New threshold for undefine cloud phase detection   
     133      real, parameter  ::  S_cld_att = 30. ! New threshold for undefine cloud phase detection
    134134
    135135! c -------------------------------------------------------
     
    144144! c -------------------------------------------------------
    145145
    146       do ic = 1, ncol
     146      DO ic = 1, ncol
    147147        pnorm_c = pnorm(:,ic,:)
    148148        where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
     
    152152        end where
    153153         x3d(:,ic,:) = x3d_c
    154 !       profSR(:,ic,:) = x3d(:,ic,:) !TIBO
    155         profSR(:,:,ic) = x3d(:,ic,:) !TIBO2
     154!    profSR(:,ic,:) = x3d(:,ic,:) !TIBO
     155    profSR(:,:,ic) = x3d(:,ic,:) !TIBO2
    156156      enddo
    157157
     
    188188      parasolrefl(:,:) = 0.0
    189189
    190       do k = 1, nrefl
    191        do ic = 1, ncol
     190      DO k = 1, nrefl
     191       DO ic = 1, ncol
    192192         parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k)
    193193       enddo
    194194      enddo
    195195
    196       do k = 1, nrefl
     196      DO k = 1, nrefl
    197197        parasolrefl(:,k) = parasolrefl(:,k) / float(ncol)
    198198! if land=1 -> parasolrefl=undef
     
    253253      srbval(5) =  7.0
    254254      srbval(6) = 10.0
    255       do i = 7, MIN(10,Nbins)
     255      DO i = 7, MIN(10,Nbins)
    256256       srbval(i) = srbval(i-1) + 5.0
    257257      enddo
     
    268268! c c- Compute CFAD
    269269! c -------------------------------------------------------
    270       do j = 1, Nlevels
    271          do ib = 1, Nbins
    272             do k = 1, Ncolumns
    273                do i = 1, Npoints
     270      DO j = 1, Nlevels
     271         DO ib = 1, Nbins
     272            DO k = 1, Ncolumns
     273               DO i = 1, Npoints
    274274                  if (x(i,k,j) /= undef) then
    275275                     if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
     
    310310      integer,parameter  ::  Ntemp=40 ! indice of the temperature vector
    311311      integer ip, k, iz, ic, ncol, nlev, i, itemp  ! loop indice
    312       real  S_cld_att ! New threshold for undefine cloud phase detection (SR=30)       
     312      real  S_cld_att ! New threshold for undefine cloud phase detection (SR=30)
    313313      integer toplvlsat  ! level of the first cloud with SR>30
    314314      real alpha50, beta50, gamma50, delta50, epsilon50, zeta50 ! Polynomial Coef of the phase
     
    316316
    317317! Input variables
    318       real tmp(Npoints,Nlevels)                 ! temperature
     318      real tmp(Npoints,Nlevels)            ! temperature
    319319      real ATB(Npoints,Ncolumns,Nlevels) ! 3D Attenuated backscatter
    320320      real ATBperp(Npoints,Ncolumns,Nlevels) ! 3D perpendicular attenuated backscatter
     
    332332
    333333! Local variables
    334       real tmpi(Npoints,Ncolumns,Nlevels)       ! temperature of ice cld
    335       real tmpl(Npoints,Ncolumns,Nlevels)       ! temperature of liquid cld
    336       real tmpu(Npoints,Ncolumns,Nlevels)       ! temperature of undef cld
     334      real tmpi(Npoints,Ncolumns,Nlevels)    ! temperature of ice cld
     335      real tmpl(Npoints,Ncolumns,Nlevels)    ! temperature of liquid cld
     336      real tmpu(Npoints,Ncolumns,Nlevels)    ! temperature of undef cld
    337337
    338338      real checktemp, ATBperp_tmp ! temporary variable
     
    399399                -54.,-51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18.,  &
    400400                -15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,21.,24.,200. /)
    401        
     401
    402402! convert C to K
    403403      tempmod=tempmod+273.15
     
    418418! ---------------------------------------------------------------
    419419
    420       do k = 1, Nlevels
     420      DO k = 1, Nlevels
    421421
    422422! cloud detection at subgrid-scale:
     
    454454      nsublay(:,:,4) = 0.0
    455455
    456       do k = Nlevels, 1, -1
    457        do ic = 1, Ncolumns
    458         do ip = 1, Npoints
     456      DO k = Nlevels, 1, -1
     457       DO ic = 1, Ncolumns
     458        DO ip = 1, Npoints
    459459
    460460         if(srok(ip,ic,k).gt.0.)then
    461461           ! Computation of the cloud fraction as a function of the temperature
    462462           ! instead of height, for ice,liquid and all clouds
    463            do itemp=1,Ntemp
     463           DO itemp=1,Ntemp
    464464             if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
    465465               lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
     
    469469
    470470         if (cldy(ip,ic,k).eq.1.) then
    471            do itemp=1,Ntemp
     471           DO itemp=1,Ntemp
    472472             if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
    473473               lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
     
    505505      cldlay = 0.0
    506506      nsublay = 0.0
    507       do k = Nlevels, 1, -1
    508        do ic = 1, Ncolumns
    509         do ip = 1, Npoints
     507      DO k = Nlevels, 1, -1
     508       DO ic = 1, Ncolumns
     509        DO ip = 1, Npoints
    510510
    511511          ! Computation of the cloud fraction as a function of the temperature
    512512          ! instead of height, for ice,liquid and all clouds
    513513          if(srok(ip,ic,k).gt.0.)then
    514           do itemp=1,Ntemp
     514          DO itemp=1,Ntemp
    515515            if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
    516516              lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
     
    520520
    521521          if(cldy(ip,ic,k).eq.1.)then
    522           do itemp=1,Ntemp
     522          DO itemp=1,Ntemp
    523523            if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
    524524              lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
     
    562562      nsublayer = 0.0
    563563
    564       do iz = 1, Ncat
    565        do ic = 1, Ncolumns
     564      DO iz = 1, Ncat
     565       DO ic = 1, Ncolumns
    566566
    567567          cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
     
    582582! 4.1 - For Cloudy pixels with 8.16km < z < 19.2km
    583583! ---------------------------------------------------------------
    584 do ncol=1,Ncolumns
    585 do i=1,Npoints
    586 
    587       do nlev=Nlevels,18,-1  ! from 19.2km until 8.16km
     584DO ncol=1,Ncolumns
     585DO i=1,Npoints
     586
     587      DO nlev=Nlevels,18,-1  ! from 19.2km until 8.16km
    588588         p1 = pplay(i,nlev)
    589589
    590590
    591591! Avoid zero values
    592         if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
     592    if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
    593593! Computation of the ATBperp along the phase discrimination line
    594594           ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
     
    609609               lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1.   ! keep the information "temperature criterium used"
    610610                                                    ! to classify the phase cloud
    611                    cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
     611                cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
    612612                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    613                    cldlayphase(i,ncol,3,2) = 1.
    614                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    615                    cldlayphase(i,ncol,2,2) = 1.
    616                 else                                                    ! low cloud
    617                    cldlayphase(i,ncol,1,2) = 1.
     613               cldlayphase(i,ncol,3,2) = 1.
     614             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     615                cldlayphase(i,ncol,2,2) = 1.
     616         else                                                    ! low cloud
     617                cldlayphase(i,ncol,1,2) = 1.
    618618                endif
    619                    cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
    620                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    621                    cldlayphase(i,ncol,3,5) = 1.
    622                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    623                    cldlayphase(i,ncol,2,5) = 1.
    624                 else                                                    ! low cloud
    625                    cldlayphase(i,ncol,1,5) = 1.
     619                cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
     620             if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     621               cldlayphase(i,ncol,3,5) = 1.
     622             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     623                cldlayphase(i,ncol,2,5) = 1.
     624         else                                                    ! low cloud
     625                cldlayphase(i,ncol,1,5) = 1.
    626626                endif
    627627
     
    630630              lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
    631631              tmpi(i,ncol,nlev)=tmp(i,nlev)
    632                    cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
    633                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    634                    cldlayphase(i,ncol,3,1) = 1.
    635                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    636                    cldlayphase(i,ncol,2,1) = 1.
    637                 else                                                    ! low cloud
    638                    cldlayphase(i,ncol,1,1) = 1.
     632                cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
     633             if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     634               cldlayphase(i,ncol,3,1) = 1.
     635             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     636                cldlayphase(i,ncol,2,1) = 1.
     637         else                                                    ! low cloud
     638                cldlayphase(i,ncol,1,1) = 1.
    639639                endif
    640640
     
    651651               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
    652652               tmpl(i,ncol,nlev)=tmp(i,nlev)
    653                    cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
    654                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    655                    cldlayphase(i,ncol,3,2) = 1. 
    656                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    657                    cldlayphase(i,ncol,2,2) = 1.
    658                 else                                                    ! low cloud
    659                    cldlayphase(i,ncol,1,2) = 1.
    660                 endif
     653                cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
     654             if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     655                cldlayphase(i,ncol,3,2) = 1.
     656             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     657                cldlayphase(i,ncol,2,2) = 1.
     658         else                                                    ! low cloud
     659                cldlayphase(i,ncol,1,2) = 1.
     660         endif
    661661
    662662             else
     
    666666               lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1.   ! keep the information "temperature criterium used"
    667667                                                    ! to classify the phase cloud
    668                    cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
    669                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    670                    cldlayphase(i,ncol,3,4) = 1. 
    671                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    672                    cldlayphase(i,ncol,2,4) = 1.
    673                 else                                                    ! low cloud
    674                    cldlayphase(i,ncol,1,4) = 1.
    675                 endif
    676                    cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
    677                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    678                    cldlayphase(i,ncol,3,1) = 1. 
    679                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    680                    cldlayphase(i,ncol,2,1) = 1.
    681                 else                                                    ! low cloud
    682                    cldlayphase(i,ncol,1,1) = 1.
    683                 endif
     668                cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
     669             if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     670                cldlayphase(i,ncol,3,4) = 1.
     671             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     672                cldlayphase(i,ncol,2,4) = 1.
     673         else                                                    ! low cloud
     674                cldlayphase(i,ncol,1,4) = 1.
     675         endif
     676                cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
     677            if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     678                cldlayphase(i,ncol,3,1) = 1.
     679             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     680                cldlayphase(i,ncol,2,1) = 1.
     681         else                                                    ! low cloud
     682                cldlayphase(i,ncol,1,1) = 1.
     683         endif
    684684
    685685             endif
    686686
    687687            endif  ! end of discrimination condition
    688         endif  ! end of cloud condition
     688    endif  ! end of cloud condition
    689689      enddo ! end of altitude loop
    690690
     
    696696
    697697      toplvlsat=0
    698       do nlev=17,1,-1  ! from 8.16km until 0km
     698      DO nlev=17,1,-1  ! from 8.16km until 0km
    699699         p1 = pplay(i,nlev)
    700700
    701         if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
     701    if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
    702702! Phase discrimination line : ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50
    703703!                                  + ATB*epsilon50 + zeta50
     
    718718               lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1.
    719719
    720                    cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
     720                cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
    721721               if ( p1.gt.0. .and. p1.lt.(440.*100.)) then              ! high cloud
    722                    cldlayphase(i,ncol,3,2) = 1.
    723                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    724                    cldlayphase(i,ncol,2,2) = 1.
    725                 else                                                    ! low cloud
    726                    cldlayphase(i,ncol,1,2) = 1.
     722               cldlayphase(i,ncol,3,2) = 1.
     723             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     724                cldlayphase(i,ncol,2,2) = 1.
     725         else                                                    ! low cloud
     726                cldlayphase(i,ncol,1,2) = 1.
    727727                endif
    728728
    729                    cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
    730                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    731                    cldlayphase(i,ncol,3,5) = 1.
    732                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    733                    cldlayphase(i,ncol,2,5) = 1.
    734                 else                                                    ! low cloud
    735                    cldlayphase(i,ncol,1,5) = 1.
     729                cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
     730             if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     731               cldlayphase(i,ncol,3,5) = 1.
     732             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     733                cldlayphase(i,ncol,2,5) = 1.
     734         else                                                    ! low cloud
     735                cldlayphase(i,ncol,1,5) = 1.
    736736                endif
    737737
     
    741741              tmpi(i,ncol,nlev)=tmp(i,nlev)
    742742
    743                    cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
    744                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    745                    cldlayphase(i,ncol,3,1) = 1.
    746                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    747                    cldlayphase(i,ncol,2,1) = 1.
    748                 else                                                    ! low cloud
    749                    cldlayphase(i,ncol,1,1) = 1.
     743                 cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
     744            if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     745               cldlayphase(i,ncol,3,1) = 1.
     746             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     747                cldlayphase(i,ncol,2,1) = 1.
     748         else                                                    ! low cloud
     749                cldlayphase(i,ncol,1,1) = 1.
    750750                endif
    751751
     
    763763               tmpl(i,ncol,nlev)=tmp(i,nlev)
    764764
    765                    cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
    766                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    767                    cldlayphase(i,ncol,3,2) = 1. 
    768                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    769                    cldlayphase(i,ncol,2,2) = 1.
    770                 else                                                    ! low cloud
    771                    cldlayphase(i,ncol,1,2) = 1.
    772                 endif
     765                cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
     766             if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     767                cldlayphase(i,ncol,3,2) = 1.
     768             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     769                cldlayphase(i,ncol,2,2) = 1.
     770         else                                                    ! low cloud
     771                cldlayphase(i,ncol,1,2) = 1.
     772         endif
    773773
    774774             else
     
    778778               lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1.  ! false liq ==> ice
    779779
    780                    cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
    781                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    782                    cldlayphase(i,ncol,3,4) = 1. 
    783                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    784                    cldlayphase(i,ncol,2,4) = 1.
    785                 else                                                    ! low cloud
    786                    cldlayphase(i,ncol,1,4) = 1.
    787                 endif
    788 
    789                    cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
    790                 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
    791                    cldlayphase(i,ncol,3,1) = 1. 
    792                 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
    793                    cldlayphase(i,ncol,2,1) = 1.
    794                 else                                                    ! low cloud
    795                    cldlayphase(i,ncol,1,1) = 1.
    796                 endif
     780                cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
     781             if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     782                cldlayphase(i,ncol,3,4) = 1.
     783             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     784                cldlayphase(i,ncol,2,4) = 1.
     785         else                                                    ! low cloud
     786                cldlayphase(i,ncol,1,4) = 1.
     787         endif
     788
     789                cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
     790            if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
     791                cldlayphase(i,ncol,3,1) = 1.
     792             else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
     793                cldlayphase(i,ncol,2,1) = 1.
     794         else                                                    ! low cloud
     795                cldlayphase(i,ncol,1,1) = 1.
     796         endif
    797797
    798798             endif
    799799           endif  ! end of discrimination condition
    800800
    801             toplvlsat=0
     801               toplvlsat=0
    802802
    803803           ! Find the level of the highest cloud with SR>30
    804             if(x(i,ncol,nlev).gt.S_cld_att)then ! SR > 30.
    805                 toplvlsat=nlev-1
    806                 goto 99
    807             endif
    808 
    809         endif  ! end of cloud condition
     804        if(x(i,ncol,nlev).gt.S_cld_att)then    ! SR > 30.
     805              toplvlsat=nlev-1
     806               goto 99
     807            endif
     808
     809    endif  ! end of cloud condition
    810810       enddo  ! end of altitude loop
    811811
     
    818818!____________________________________________________________________________________________________
    819819
    820 if(toplvlsat.ne.0)then         
    821       do nlev=toplvlsat,1,-1
     820if(toplvlsat.ne.0)then
     821      DO nlev=toplvlsat,1,-1
    822822         p1 = pplay(i,nlev)
    823         if(cldy(i,ncol,nlev).eq.1.)then
     823    if(cldy(i,ncol,nlev).eq.1.)then
    824824           lidarcldphase(i,nlev,3)=lidarcldphase(i,nlev,3)+1.
    825825           tmpu(i,ncol,nlev)=tmp(i,nlev)
    826826
    827                    cldlayphase(i,ncol,4,3) = 1.                         ! tot cloud
     827                cldlayphase(i,ncol,4,3) = 1.                         ! tot cloud
    828828          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then              ! high cloud
    829829             cldlayphase(i,ncol,3,3) = 1.
    830830          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid cloud
    831831             cldlayphase(i,ncol,2,3) = 1.
    832           else                                                     ! low cloud
     832      else                                                     ! low cloud
    833833             cldlayphase(i,ncol,1,3) = 1.
    834           endif
    835 
    836         endif   
     834      endif
     835
     836        endif
    837837      enddo
    838838endif
     
    876876
    877877! Compute Phase low mid high cloud fractions
    878     do iz = 1, Ncat
    879        do i=1,Nphase-3
    880        do ic = 1, Ncolumns
     878    DO iz = 1, Ncat
     879       DO i=1,Nphase-3
     880       DO ic = 1, Ncolumns
    881881          cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
    882882          cldlayerphasesum(:,iz)=cldlayerphasesum(:,iz)+cldlayphase(:,ic,iz,i)
     
    885885    enddo
    886886
    887     do iz = 1, Ncat
    888        do i=4,5
    889        do ic = 1, Ncolumns
     887    DO iz = 1, Ncat
     888       DO i=4,5
     889       DO ic = 1, Ncolumns
    890890          cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)         
    891891       enddo
     
    901901    ENDWHERE
    902902
    903     do i=1,Nphase-1
     903    DO i=1,Nphase-1
    904904      WHERE ( cldlayerphasesum(:,:).gt.0.0 )
    905905         cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
     
    908908
    909909
    910     do i=1,Npoints
    911        do iz=1,Ncat
     910    DO i=1,Npoints
     911       DO iz=1,Ncat
    912912          checkcldlayerphase=0.
    913913          checkcldlayerphase2=0.
    914914
    915915          if (cldlayerphasesum(i,iz).gt.0.0 )then
    916              do ic=1,Nphase-3
     916             DO ic=1,Nphase-3
    917917                checkcldlayerphase=checkcldlayerphase+cldlayerphase(i,iz,ic) 
    918918             enddo
     
    925925    enddo
    926926
    927     do i=1,Nphase-1
     927    DO i=1,Nphase-1
    928928      WHERE ( nsublayer(:,:).eq.0.0 )
    929929         cldlayerphase(:,:,i) = undef
     
    934934
    935935! Compute Phase 3D as a function of temperature
    936 do nlev=1,Nlevels
    937 do ncol=1,Ncolumns     
    938 do i=1,Npoints
    939 do itemp=1,Ntemp
     936DO nlev=1,Nlevels
     937DO ncol=1,Ncolumns
     938DO i=1,Npoints
     939DO itemp=1,Ntemp
    940940if(tmpi(i,ncol,nlev).gt.0.)then
    941941      if( (tmpi(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpi(i,ncol,nlev).lt.tempmod(itemp+1)) )then
     
    957957
    958958! Check temperature cloud fraction
    959 do i=1,Npoints
    960    do itemp=1,Ntemp
     959DO i=1,Npoints
     960   DO itemp=1,Ntemp
    961961checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4)
    962962
    963         if(checktemp.NE.lidarcldtemp(i,itemp,1))then
    964           print *, i,itemp
    965           print *, lidarcldtemp(i,itemp,1:4)
    966         endif
     963    if(checktemp.NE.lidarcldtemp(i,itemp,1))then
     964      print *, i,itemp
     965      print *, lidarcldtemp(i,itemp,1:4)
     966    endif
    967967
    968968   enddo
     
    979979ENDWHERE
    980980
    981 do i=1,4
     981DO i=1,4
    982982  WHERE(lidarcldtempind(:,:).gt.0.)
    983983     lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
     
    10251025 
    10261026    ! ####################################################################################
    1027         ! 1) Initialize   
     1027    ! 1) Initialize
    10281028    ! ####################################################################################
    10291029    cldtype               = 0.0
     
    10401040    ! 2) Cloud detection and Fully attenuated layer detection
    10411041    ! ####################################################################################
    1042     do k=1,Nlevels
     1042    DO k=1,Nlevels
    10431043       ! Cloud detection at subgrid-scale:
    10441044       where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
     
    10721072    ! ####################################################################################
    10731073
    1074     do k= Nlevels,1,-1
    1075        do ic = 1, Ncolumns
    1076           do ip = 1, Npoints
     1074    DO k= Nlevels,1,-1
     1075       DO ic = 1, Ncolumns
     1076          DO ip = 1, Npoints
    10771077
    10781078             cldlay(ip,ic,1)   = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque clouds
     
    10901090
    10911091! OPAQ variables
    1092      do ic = 1, Ncolumns
    1093         do ip = 1, Npoints
     1092     DO ic = 1, Ncolumns
     1093        DO ip = 1, Npoints
    10941094
    10951095     ! Declaring non-opaque cloudy profiles as thin cloud profiles
    1096            if ( (cldlay(ip,ic,4) .eq. 1.0) .and. (cldlay(ip,ic,1) .eq. 0.0) ) then
    1097               cldlay(ip,ic,2)  =  1.0
    1098            endif
     1096       if ( (cldlay(ip,ic,4) .eq. 1.0) .and. (cldlay(ip,ic,1) .eq. 0.0) ) then
     1097          cldlay(ip,ic,2)  =  1.0
     1098        endif
    10991099
    11001100     ! Filling in 3D and 2D variables
    11011101
    11021102     ! Opaque cloud profiles
    1103            if ( cldlay(ip,ic,1) .eq. 1.0 ) then
    1104               zopac = 0.0
    1105               do k=2,Nlevels
     1103       if ( cldlay(ip,ic,1) .eq. 1.0 ) then
     1104          zopac = 0.0
     1105          DO k=2,Nlevels
    11061106     ! Declaring opaque cloud fraction and z_opaque altitude for 3D and 2D variables
    1107                  if ( (cldy(ip,ic,k) .eq. 1.0) .and. (zopac .eq. 0.0) ) then
    1108                     lidarcldtype(ip,k-1,3) = lidarcldtype(ip,k-1,3) + 1.0
    1109                     cldlay(ip,ic,3)        = vgrid_z(k-1) !z_opaque altitude
    1110                     nsublay(ip,ic,3)       = 1.0
    1111                     zopac = 1.0
    1112                 endif
    1113                  if ( cldy(ip,ic,k) .eq. 1.0 ) then
    1114                     lidarcldtype(ip,k,1)   = lidarcldtype(ip,k,1) + 1.0
     1107             if ( (cldy(ip,ic,k) .eq. 1.0) .and. (zopac .eq. 0.0) ) then
     1108            lidarcldtype(ip,k-1,3) = lidarcldtype(ip,k-1,3) + 1.0
     1109            cldlay(ip,ic,3)        = vgrid_z(k-1) !z_opaque altitude
     1110            nsublay(ip,ic,3)       = 1.0
     1111            zopac = 1.0
     1112        endif
     1113             if ( cldy(ip,ic,k) .eq. 1.0 ) then
     1114            lidarcldtype(ip,k,1)   = lidarcldtype(ip,k,1) + 1.0
    11151115                 endif
    1116               enddo
    1117            endif
     1116          enddo
     1117       endif
    11181118
    11191119     ! Thin cloud profiles
    1120            if ( cldlay(ip,ic,2) .eq. 1.0 ) then
    1121               do k=1,Nlevels
     1120       if ( cldlay(ip,ic,2) .eq. 1.0 ) then
     1121          DO k=1,Nlevels
    11221122     ! Declaring thin cloud fraction for 3D variable
    11231123                 if ( cldy(ip,ic,k) .eq. 1.0 ) then
    11241124                    lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1.0
    11251125                 endif
    1126               enddo
     1126          enddo
    11271127           endif
    11281128
     
    11461146    ! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=Nlevels) to SFC(k=1)
    11471147       lidarcldtype(:,Nlevels,4) = lidarcldtype(:,Nlevels,3)
    1148     do ip = 1, Npoints
    1149         do k = Nlevels-1, 1, -1
     1148    DO ip = 1, Npoints
     1149         DO k = Nlevels-1, 1, -1
    11501150           if ( lidarcldtype(ip,k,3) .ne. undef ) then
    1151               lidarcldtype(ip,k,4) = lidarcldtype(ip,k+1,4) + lidarcldtype(ip,k,3)
     1151          lidarcldtype(ip,k,4) = lidarcldtype(ip,k+1,4) + lidarcldtype(ip,k,3)
    11521152           endif
    1153         enddo
     1153    enddo
    11541154    enddo
    11551155    where ( nsubopaq(:,:) .eq. 0.0 )
     
    11591159    ! Layered cloud types (opaque, thin and z_opaque 2D variables)
    11601160
    1161     do iz = 1, Ntype
    1162        do ic = 1, Ncolumns
     1161    DO iz = 1, Ntype
     1162       DO ic = 1, Ncolumns
    11631163          cldtype(:,iz)   = cldtype(:,iz)   + cldlay(:,ic,iz)
    11641164          nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_modis_sim.F90

    r5099 r5158  
    200200    logical, dimension(size(retrievedTau))                     :: cloudMask
    201201    real,    dimension(size(waterSize, 1), size(waterSize, 2)) :: tauLiquidFraction, tauTotal
    202     real    :: integratedLiquidFraction
     202    REAL    :: integratedLiquidFraction
    203203    integer :: i, nSubcols, nLevels
    204204
     
    256256    end if
    257257   
    258     do i = 1, nSubCols
     258    DO i = 1, nSubCols
    259259      if(cloudMask(i)) then
    260260
     
    371371           liquid_opticalThickness, ice_opticalThickness, tauLiquidFraction
    372372   
    373     real :: seuil
     373    REAL :: seuil
    374374    ! ---------------------------------------------------
    375375   
     
    596596    reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
    597597    reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
    598     do j=1,nPoints
     598    DO j=1,nPoints
    599599
    600600       ! Fill clear and optically thin subcolumns with fill
     
    651651    integer :: ij,ik
    652652   
    653     do ij=2,nbin1+1
    654        do ik=2,nbin2+1
     653    DO ij=2,nbin1+1
     654       DO ik=2,nbin2+1
    655655          jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
    656656               var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))       
     
    778778    ! Joint histogram
    779779
    780     do i = 1, numTauHistogramBins
     780    DO i = 1, numTauHistogramBins
    781781      where(cloudMask(:, :))
    782782        tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. &
     
    787787    end do
    788788
    789     do i = 1, numPressureHistogramBins
     789    DO i = 1, numPressureHistogramBins
    790790      where(cloudMask(:, :))
    791791        pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. &
     
    796796    end do
    797797   
    798     do i = 1, numPressureHistogramBins
    799       do j = 1, numTauHistogramBins
     798    DO i = 1, numPressureHistogramBins
     799      DO j = 1, numTauHistogramBins
    800800        Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = &
    801801          real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols)
     
    808808    real, dimension(:), intent(in) :: tauIncrement, pressure
    809809    real,               intent(in) :: tauLimit
    810     real                           :: cloud_top_pressure
     810    REAL                           :: cloud_top_pressure
    811811
    812812    ! Find the extinction-weighted pressure. Assume that pressure varies linearly between
    813813    !   layers and use the trapezoidal rule.
    814814
    815     real :: deltaX, totalTau, totalProduct
     815    REAL :: deltaX, totalTau, totalProduct
    816816    integer :: i
    817817   
    818818    totalTau = 0.; totalProduct = 0.
    819     do i = 2, size(tauIncrement)
     819    DO i = 2, size(tauIncrement)
    820820      if(totalTau + tauIncrement(i) > tauLimit) then
    821821        deltaX = tauLimit - totalTau
     
    840840    real, dimension(:), intent(in) :: tauIncrement, f
    841841    real,               intent(in) :: tauLimit
    842     real                           :: weight_by_extinction
     842    REAL                           :: weight_by_extinction
    843843
    844844    ! Find the extinction-weighted value of f(tau), assuming constant f within each layer
    845845
    846     real    :: deltaX, totalTau, totalProduct
     846    REAL    :: deltaX, totalTau, totalProduct
    847847    integer :: i
    848848   
    849849    totalTau = 0.; totalProduct = 0.
    850     do i = 1, size(tauIncrement)
     850    DO i = 1, size(tauIncrement)
    851851      if(totalTau + tauIncrement(i) > tauLimit) then
    852852        deltaX       = tauLimit - totalTau
     
    864864  pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size)
    865865    real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size
    866     real                           :: compute_nir_reflectance
     866    REAL                           :: compute_nir_reflectance
    867867   
    868868    real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, &
     
    893893      integer, intent(in) :: phase
    894894      real,    intent(in) :: tau, obs_Refl_nir
    895       real                :: retrieve_re
     895      REAL                :: retrieve_re
    896896
    897897      ! Finds the re that produces the minimum mis-match between predicted and observed reflectance in
     
    907907
    908908      real, parameter :: min_distance_to_boundary = 0.01
    909       real    :: re_min, re_max, delta_re
     909      REAL    :: re_min, re_max, delta_re
    910910      integer :: i
    911911     
     
    957957    real, dimension(:), intent(in) :: x, y
    958958    real,               intent(in) :: yobs
    959     real                           :: interpolate_to_min
     959    REAL                           :: interpolate_to_min
    960960
    961961    ! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs)
     
    10081008    integer, intent(in) :: phase
    10091009    real,    intent(in) :: re
    1010     real :: get_g_nir
     1010    REAL :: get_g_nir
    10111011
    10121012    real, dimension(3), parameter :: ice_coefficients         = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
     
    10351035        integer, intent(in) :: phase
    10361036        real,    intent(in) :: re
    1037         real                :: get_ssa_nir
     1037        REAL                :: get_ssa_nir
    10381038
    10391039        ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
     
    10601060    real,               intent(in) :: x
    10611061    real, dimension(:), intent(in) :: coefficients
    1062     real                           :: fit_to_cubic
     1062    REAL                           :: fit_to_cubic
    10631063   
    10641064   
     
    10691069    real,               intent(in) :: x
    10701070    real, dimension(:), intent(in) :: coefficients
    1071     real                           :: fit_to_quadratic
     1071    REAL                           :: fit_to_quadratic
    10721072   
    10731073   
     
    10791079  pure function compute_toa_reflectace(tau, g, w0)
    10801080    real, dimension(:), intent(in) :: tau, g, w0
    1081     real                           :: compute_toa_reflectace
     1081    REAL                           :: compute_toa_reflectace
    10821082   
    10831083    logical, dimension(size(tau))         :: cloudMask
    10841084    integer, dimension(count(tau(:) > 0)) :: cloudIndicies
    10851085    real,    dimension(count(tau(:) > 0)) :: Refl,     Trans
    1086     real                                  :: Refl_tot, Trans_tot
     1086    REAL                                  :: Refl_tot, Trans_tot
    10871087    integer                               :: i
    10881088    ! ---------------------------------------
     
    10921092    cloudMask = tau(:) > 0.
    10931093    cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask)
    1094     do i = 1, size(cloudIndicies)
     1094    DO i = 1, size(cloudIndicies)
    10951095      call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
    10961096    end do
     
    11151115    integer, parameter :: beam = 2
    11161116    real,    parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
    1117     real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
     1117    REAL :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
    11181118            rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th
    11191119
     
    11831183  elemental function two_stream_reflectance(tauint, gint, w0int)
    11841184    real, intent(in) :: tauint, gint, w0int
    1185     real             :: two_stream_reflectance
     1185    REAL             :: two_stream_reflectance
    11861186
    11871187    ! Compute reflectance in a single layer using the two stream approximation
     
    11941194    integer, parameter :: beam = 2
    11951195    real,    parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
    1196     real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
     1196    REAL :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
    11971197            rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
    11981198    ! ------------------------
     
    12671267      Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1)   
    12681268     
    1269       do i=2, size(Refl)
     1269      DO i=2, size(Refl)
    12701270          ! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
    12711271          Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1 - Refl_cumulative(i-1) * Refl(i))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/optics_lib.F90

    r5099 r5158  
    519519
    520520!   // region from 0.045 microns to 167.0 microns - no temperature depend
    521     do i=2,nwl
     521    DO i=2,nwl
    522522      if(alam < wl(i)) continue
    523523    enddo
     
    539539    if(tk > temref(1)) tk=temref(1)
    540540    if(tk < temref(4)) tk=temref(4)
    541     do 11 i=2,4
     541    DO 11 i=2,4
    542542      if(tk.ge.temref(i)) go to 12
    543543    11 continue
    544544    12 lt1=i
    545545    lt2=i-1
    546     do 13 i=2,nwlt
     546    DO 13 i=2,nwlt
    547547      if(alam.le.wlt(i)) go to 14
    548548    13 continue
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/phys_cosp.F90

    r5133 r5158  
    4444! mr_ozone,                             !Concentration ozone (Kg/Kg)
    4545! dem_s                                 !Cloud optical emissivity
    46 ! dtau_s                                !Cloud optical thickness
    47 ! emsfc_lw = 1.                         !Surface emissivity dans radlwsw.F90
     46! dtau_s                           !Cloud optical thickness
     47! emsfc_lw = 1.                    !Surface emissivity dans radlwsw.F90
    4848
    4949!!! Outputs :
     
    132132! Declaration necessaires pour les sorties IOIPSL
    133133  integer :: ii
    134   real    :: ecrit_day,ecrit_hf,ecrit_mth, missing_val
     134  REAL    :: ecrit_day,ecrit_hf,ecrit_mth, missing_val
    135135  logical :: ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, ok_all_xml
    136136
     
    150150  real,dimension(Nlevlmdz)        :: presnivs
    151151  integer                         :: itap,k,ip
    152   real                            :: dtime,freq_cosp
     152  REAL                            :: dtime,freq_cosp
    153153  real,dimension(2)               :: time_bnds
    154154
     
    249249
    250250        zlev_half(:,1) = phis(:)/9.81
    251         do k = 2, Nlevels
    252           do ip = 1, Npoints
     251        DO k = 2, Nlevels
     252          DO ip = 1, Npoints
    253253           zlev_half(ip,k) = phi(ip,k)/9.81 + &
    254254               (phi(ip,k)-phi(ip,k-1))/9.81 * (ph(ip,k)-p(ip,k)) / (p(ip,k)-p(ip,k-1))
     
    267267        gbx%skt  = skt !Skin temperature (K)
    268268
    269         do ip = 1, Npoints
     269        DO ip = 1, Npoints
    270270          if (fracTerLic(ip).ge.0.5) then
    271271             gbx%land(ip) = 1.
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator.F90

    r5099 r5158  
    4141!   [hm_matrix]       table of hydrometeor mixing rations (g/kg)
    4242!   [re_matrix]       table of hydrometeor effective radii.  0 ==> use defaults. (units=microns)   
    43 !   [Np_matrix]       table of hydrometeor number concentration.  0 ==> use defaults. (units = 1/kg)
     43!   [Np_matrix]          table of hydrometeor number concentration.  0 ==> use defaults. (units = 1/kg)
    4444
    4545! Outputs:
     
    117117  ns       ! number of discrete drop sizes
    118118
    119   logical :: hydro      ! true=hydrometeor in vol, false=none
     119  logical :: hydro    ! true=hydrometeor in vol, false=none
    120120  real*8 :: &
    121121  rho_a, &   ! air density (kg m^-3)
     
    180180    d_gate=-1
    181181  endif
    182   do k=start_gate,end_gate,d_gate
     182  DO k=start_gate,end_gate,d_gate
    183183  ! // loop over each profile (nprof)
    184     do pr=1,nprof
     184    DO pr=1,nprof
    185185      t_kelvin = t_matrix(pr,k)
    186186!     :: determine if hydrometeor(s) present in volume
    187187      hydro = .false.
    188       do j=1,hp%nhclass
     188      DO j=1,hp%nhclass
    189189        if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then
    190190          hydro = .true.
     
    198198        rho_a = (p_matrix(pr,k)*100.)/(287.0*(t_kelvin))
    199199!       :: loop over hydrometeor type
    200         do tp=1,hp%nhclass
     200        DO tp=1,hp%nhclass
    201201          if (hm_matrix(tp,pr,k) <= 1E-12) cycle
    202202          phase = hp%phase(tp)
     
    233233            if (iRe_type.lt.1) iRe_type=1
    234234
    235             Re=step*(iRe_type+0.5)      ! set value of Re to closest value allowed in LUT.
     235            Re=step*(iRe_type+0.5)    ! set value of Re to closest value allowed in LUT.
    236236            iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step)
    237237
     
    273273              if (hp%rho(tp) < 0) then
    274274                ! Use equivalent volume spheres.
    275                 hp%rho_eff(tp,1:ns,iRe_type) = 917                              ! solid ice == equivalent volume approach
     275                hp%rho_eff(tp,1:ns,iRe_type) = 917                  ! solid ice == equivalent volume approach
    276276                Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * ( (Di*1E-6) ** (hp%bpm(tp)/3.0) )  * 1E6
    277277                ! alternative is to comment out above two lines and use the following block
     
    362362          endif
    363363
    364         enddo   ! end loop of tp (hydrometeor type)
     364        enddo    ! end loop of tp (hydrometeor type)
    365365
    366366      else
    367 !     :: volume is hydrometeor-free     
     367!     :: volume is hydrometeor-free
    368368        kr_vol(pr,k) = 0
    369369        z_vol(pr,k)  = undef
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator_init.F90

    r5099 r5158  
    6060   integer, intent(in) :: nhclass       ! number of hydrometeor classes in use
    6161   integer, intent(in) :: use_gas_abs,do_ray
    62    real :: undef
     62   REAL :: undef
    6363   real,dimension(nhclass),intent(in)     ::    hclass_dmin,hclass_dmax, &
    6464                                hclass_apm,hclass_bpm,hclass_rho, &
     
    9090        ! Store settings for hydrometeor types in hp (class_parameter) structure.   
    9191
    92         do i = 1,nhclass
     92        DO i = 1,nhclass
    9393        hp%dtype(i) = hclass_type(i)
    9494        hp%phase(i) = hclass_phase(i)
     
    117117
    118118    hp%base_list(1)=0;
    119     do j=1,Re_MAX_BIN
     119    DO j=1,Re_MAX_BIN
    120120        hp%step_list(j)=0.1+0.1*((j-1)**1.5);
    121121        if(hp%step_list(j)>Re_BIN_LENGTH) then
     
    129129    ! set up Temperature bin structure used for z scaling
    130130
    131     do i=1,cnt_ice
     131    DO i=1,cnt_ice
    132132        hp%mt_tti(i)=(i-1)*5-90 + 273.15
    133133    enddo
    134134   
    135     do i=1,cnt_liq
     135    DO i=1,cnt_liq
    136136        hp%mt_ttl(i)=(i-1)*5-60 + 273.15
    137137    enddo
     
    143143
    144144        hp%D(1) = dmin
    145         do i=2,nd
     145        DO i=2,nd
    146146          hp%D(i) = hp%D(i-1)*deltp
    147147        enddo   
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator_types.F90

    r5095 r5158  
    4545   
    4646    ! defines location of radar relative to hgt_matrix.   
    47     logical :: radar_at_layer_one       ! if true radar is assume to be at the edge
    48                                         ! of the first layer, if the first layer is the
    49                                         ! surface than a ground-based radar.   If the
    50                                         ! first layer is the top-of-atmosphere, then
    51                                         ! a space borne radar.
     47    logical :: radar_at_layer_one    ! if true radar is assume to be at the edge
     48                        ! of the first layer, if the first layer is the
     49                        ! surface than a ground-based radar.   If the
     50                        ! first layer is the top-of-atmosphere, then
     51                        ! a space borne radar.
    5252   
    5353    ! variables used to store Z scale factors
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/scale_luts_io.F90

    r5099 r5158  
    4242        trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
    4343   
    44             do i=1,maxhclass
    45             do j=1,mt_ntt
    46             do k=1,nRe_types
     44            DO i=1,maxhclass
     45            DO j=1,mt_ntt
     46            DO k=1,nRe_types
    4747       
    4848            ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt)
     
    9595        trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat'
    9696   
    97         do i=1,maxhclass
    98         do j=1,mt_ntt
    99         do k=1,nRe_types
     97        DO i=1,maxhclass
     98        DO j=1,mt_ntt
     99        DO k=1,nRe_types
    100100       
    101101            ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/zeff.F90

    r5099 r5158  
    110110    sizep = (pi*D0)/wl
    111111    dqv(1) = 0.
    112     do i=1,nsizes
     112    DO i=1,nsizes
    113113      call mieint(sizep(i), m0(i), one, dqv, qext(i), dqsc, qbsca(i), &
    114114        dg, xs1, xs2, dph, err)
Note: See TracChangeset for help on using the changeset viewer.