Ignore:
Timestamp:
Aug 2, 2024, 2:12:03 PM (6 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/ecrad
Files:
55 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver.F90

    r4853 r5158  
    311311  tstart = omp_get_wtime()
    312312#endif
    313   do jrepeat = 1,driver_config%nrepeat
     313  DO jrepeat = 1,driver_config%nrepeat
    314314   
    315315    if (driver_config%do_parallel) then
     
    321321     
    322322      !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME)
    323       do jblock = 1, nblock
     323      DO jblock = 1, nblock
    324324        ! Specify the range of columns to process.
    325325        istartcol = (jblock-1) * driver_config%nblocksize &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_config.F90

    r4773 r5158  
    400400    ! wavelength bound, noting that the number of diagnostics is one
    401401    ! fewer than the number of valid bounds
    402     do jdiag = 0,NMaxSpectralDiag
     402    DO jdiag = 0,NMaxSpectralDiag
    403403      if (this%sw_diag_wavelength_bound(jdiag+1) < 0.0_jprb) then
    404404        this%n_sw_diag = max(0,jdiag-1)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_read_input.F90

    r4773 r5158  
    563563
    564564    ! Loop through all radiatively important gases
    565     do jgas = 1,NMaxGases
     565    DO jgas = 1,NMaxGases
    566566      if (jgas == IH2O) then
    567567        if (file%exists('q')) then
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver.F90

    r4773 r5158  
    369369  tstart = omp_get_wtime()
    370370#endif
    371   do jrepeat = 1,driver_config%nrepeat
     371  DO jrepeat = 1,driver_config%nrepeat
    372372
    373373!    if (driver_config%do_parallel) then
     
    379379
    380380      !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME)
    381       do jblock = 1, nblock
     381      DO jblock = 1, nblock
    382382        ! Specify the range of columns to process.
    383383        istartcol = (jblock-1) * driver_config%nblocksize &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver_blocked.F90

    r4773 r5158  
    371371  tstart = omp_get_wtime()
    372372#endif
    373   do jrepeat = 1,driver_config%nrepeat
     373  DO jrepeat = 1,driver_config%nrepeat
    374374
    375375!    if (driver_config%do_parallel) then
     
    378378      !$OMP PARALLEL DO SCHEDULE(DYNAMIC,1)&
    379379      !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB)
    380       do jrl=1,ncol,nproma
     380      DO jrl=1,ncol,nproma
    381381        ibeg=jrl
    382382        iend=min(ibeg+nproma-1,ncol)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ifs_blocking.F90

    r4773 r5158  
    329329  !$OMP PARALLEL DO SCHEDULE(RUNTIME)&
    330330  !$OMP&PRIVATE(IB,IFLD)
    331   do ib=1,ngpblks
    332     do ifld=1,ifs_config%ifldstot
     331  DO ib=1,ngpblks
     332    DO ifld=1,ifs_config%ifldstot
    333333      zrgp(:,ifld,ib) = 0._jprb
    334334    enddo
     
    351351    !$OMP PARALLEL DO SCHEDULE(RUNTIME)&
    352352    !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB,JAER,JOFF,JLEV,JALB)
    353     do jrl=1,ncol,nproma
     353    DO jrl=1,ncol,nproma
    354354
    355355      ibeg=jrl
     
    363363      zrgp(1:il,ifs_config%iamu0,ib)  =  single_level%cos_sza(ibeg:iend)   ! cosine of solar zenith ang (mu0)
    364364
    365       do jemiss=1,yderad%nlwemiss
     365      DO jemiss=1,yderad%nlwemiss
    366366        zrgp(1:il,ifs_config%iemiss+jemiss-1,ib)  =  single_level%lw_emissivity(ibeg:iend,jemiss)
    367367      enddo
     
    378378      ! zrgp(1:il,islon,ib)    = ???
    379379
    380       do jalb=1,yderad%nsw
     380      DO jalb=1,yderad%nsw
    381381        zrgp(1:il,ifs_config%iald+jalb-1,ib)  =  single_level%sw_albedo(ibeg:iend,jalb)
    382382      enddo
    383383
    384384      if (allocated(single_level%sw_albedo_direct)) then
    385         do jalb=1,yderad%nsw
     385        DO jalb=1,yderad%nsw
    386386          zrgp(1:il,ifs_config%ialp+jalb-1,ib)  =  single_level%sw_albedo_direct(ibeg:iend,jalb)
    387387        end do
    388388      else
    389         do jalb=1,yderad%nsw
     389        DO jalb=1,yderad%nsw
    390390          zrgp(1:il,ifs_config%ialp+jalb-1,ib)  =  single_level%sw_albedo(ibeg:iend,jalb)
    391391        end do
    392392      end if
    393393     
    394       do jlev=1,nlev
     394      DO jlev=1,nlev
    395395        zrgp(1:il,ifs_config%iti+jlev-1,ib)   = temperature_fl(ibeg:iend,jlev) ! full level temperature
    396396        zrgp(1:il,ifs_config%ipr+jlev-1,ib)   = pressure_fl(ibeg:iend,jlev) ! full level pressure
     
    398398      enddo
    399399
    400       do jlev=1,nlev
     400      DO jlev=1,nlev
    401401        zrgp(1:il,ifs_config%iwv+jlev-1,ib)   = gas%mixing_ratio(ibeg:iend,jlev,IH2O) ! this is already in MassMixingRatio units
    402402        if (rad_config%do_clouds) then
     
    421421      if (yderad%naermacc == 1) then
    422422        joff=ifs_config%iaero
    423         do jaer=1,rad_config%n_aerosol_types
    424           do jlev=1,nlev
     423        DO jaer=1,rad_config%n_aerosol_types
     424          DO jlev=1,nlev
    425425            zrgp(1:il,joff,ib) = aerosol%mixing_ratio(ibeg:iend,jlev,jaer)
    426426            joff=joff+1
     
    429429      endif
    430430
    431       do jlev=1,nlev+1
     431      DO jlev=1,nlev+1
    432432        ! zrgp(1:il,ihpr+jlev-1,ib)  = ???
    433433        zrgp(1:il,ifs_config%iaprs+jlev-1,ib) = thermodynamics%pressure_hl(ibeg:iend,jlev)
     
    453453      ! local workaround variables for standalone input files
    454454      if (rad_config%do_clouds) then
    455         do jlev=1,nlev
     455        DO jlev=1,nlev
    456456          ! missing full-level temperature and pressure as well as land-sea-mask
    457457          zrgp(1:il,ifs_config%ire_liq+jlev-1,ib) = cloud%re_liq(ibeg:iend,jlev)
    458458          zrgp(1:il,ifs_config%ire_ice+jlev-1,ib) = cloud%re_ice(ibeg:iend,jlev)
    459459        enddo
    460         do jlev=1,nlev-1
     460        DO jlev=1,nlev-1
    461461          ! for the love of it, I can't figure this one out. Probably to do with
    462462          ! my crude approach of setting PGEMU?
     
    465465        if(present(iseed)) iseed(1:il,ib) = single_level%iseed(ibeg:iend)
    466466      else
    467         do jlev=1,nlev
     467        DO jlev=1,nlev
    468468          ! missing full-level temperature and pressure as well as land-sea-mask
    469469          zrgp(1:il,ifs_config%ire_liq+jlev-1,ib) = 0._jprb
    470470          zrgp(1:il,ifs_config%ire_ice+jlev-1,ib) = 0._jprb
    471471        enddo
    472         do jlev=1,nlev-1
     472        DO jlev=1,nlev-1
    473473          zrgp(1:il,ifs_config%ioverlap+jlev-1,ib) = 0._jprb
    474474        enddo
     
    531531    !$OMP PARALLEL DO SCHEDULE(RUNTIME)&
    532532    !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB,JLEV,JG)
    533     do jrl=1,ncol,nproma
     533    DO jrl=1,ncol,nproma
    534534      ibeg=jrl
    535535      iend=min(ibeg+nproma-1,ncol)
     
    537537      ib=(jrl-1)/nproma+1
    538538
    539       do jlev=1,nlev+1
     539      DO jlev=1,nlev+1
    540540        flux%sw_up(ibeg:iend,jlev) = zrgp(1:il,ifs_config%ifrso+jlev-1,ib)
    541541        flux%lw_up(ibeg:iend,jlev) = zrgp(1:il,ifs_config%ifrth+jlev-1,ib)
     
    561561      emissivity_out(ibeg:iend) = zrgp(1:il,ifs_config%iemit,ib)
    562562      if (yradiation%yrerad%lapproxswupdate) then
    563         do jg=1,yradiation%yrerad%nsw
     563        DO jg=1,yradiation%yrerad%nsw
    564564          flux_diffuse_band(ibeg:iend,jg) = zrgp(1:il,ifs_config%iswdiffuseband+jg-1,ib)
    565565          flux_direct_band(ibeg:iend,jg) = zrgp(1:il,ifs_config%iswdirectband+jg-1,ib)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_cloud_generator.F90

    r4773 r5158  
    5757       &  od_scaling, total_cloud_cover)
    5858
    59   do jlev = 1,nlev
    60     do jcol = 1,ncol
     59  DO jlev = 1,nlev
     60    DO jcol = 1,ncol
    6161      write(6,'(f5.2,a)',advance='no') od_scaling(jcol,jlev), ' '
    6262    end do
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_random_number_generator.F90

    r4773 r5158  
    2929  call random_number_generator%initialize(itype=IRngType, iseed=212075152, &
    3030       &                                  nmaxstreams=streammax)
    31   do jl = 1,1
     31  DO jl = 1,1
    3232    print *, 'initial_state = [ ', int(random_number_generator%istate(1:streammax),jpib), ' ]'
    3333    call random_number_generator%uniform_distribution(vec)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_spartacus_math.F90

    r4773 r5158  
    4545
    4646  write(*,*) 'A ='
    47   do j = 1,m
     47  DO j = 1,m
    4848     write(*,*) A(1,j,:)
    4949  end do
    5050
    5151  write(*,*) 'B ='
    52   do j = 1,m
     52  DO j = 1,m
    5353     write(*,*) B(1,j,:)
    5454  end do
     
    6060
    6161  write(*,*) 'C = A*B ='
    62   do j = 1,m
     62  DO j = 1,m
    6363     write(*,*) C(1,j,:)
    6464  end do
     
    7070
    7171  write(*,*) 'C = A\B ='
    72   do j = 1,m
     72  DO j = 1,m
    7373     write(*,*) C(1,j,:)
    7474  end do
     
    7878
    7979  write(*,*) 'expm(A) ='
    80   do j = 1,m
     80  DO j = 1,m
    8181     write(*,*) A(1,j,:)
    8282  end do
     
    9393   
    9494    write(*,*) 'fast_expm(A) = '
    95     do j = 1,m
     95    DO j = 1,m
    9696      write(*,*) A(1,j,:)
    9797    end do
     
    106106   
    107107    write(*,*) 'expm(A) = '
    108     do j = 1,m
     108    DO j = 1,m
    109109      write(*,*) A(1,j,:)
    110110    end do
     
    115115   
    116116    write(*,*) 'expm(zeros) = '
    117     do j = 1,m
     117    DO j = 1,m
    118118      write(*,*) A(1,j,:)
    119119    end do
     
    124124   
    125125    write(*,*) 'fast_expm(A) = '
    126     do j = 1,m
     126    DO j = 1,m
    127127      write(*,*) A(1,j,:)
    128128    end do
     
    140140   
    141141    write(*,*) 'expm(A) = '
    142     do j = 1,m
     142    DO j = 1,m
    143143      write(*,*) A(1,j,:)
    144144    end do
     
    149149   
    150150    write(*,*) 'expm(zeros) = '
    151     do j = 1,m
     151    DO j = 1,m
    152152      write(*,*) A(1,j,:)
    153153    end do
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/ifs/fcttre.func.h

    r4876 r5158  
    170170FOEELIQ2ICE( PTARE )  = EXP(      (PTARE-RTT)*(R3LES/(PTARE-R4LES) - R3IES/(PTARE-R4IES)))
    171171FOELSON( PTARE ) = EXP( -6096.9385_JPRB/PTARE + 21.2409642_JPRB &
    172                      - 2.711193E-2_JPRB * PTARE    &
     172                 - 2.711193E-2_JPRB * PTARE    &
    173173                     + 1.673952E-5_JPRB * PTARE**2 &
    174                      + 2.433502_JPRB * LOG(PTARE))
     174             + 2.433502_JPRB * LOG(PTARE))
    175175
    176176REAL(KIND=JPRB) :: FOEEWM_V,FOEEWMCU_V,FOELES_V,FOEIES_V
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.F90

    r4853 r5158  
    5454
    5555      ! Loop over level
    56      do jlev = istartlev,iendlev
     56     DO jlev = istartlev,iendlev
    5757
    5858        ! Loop over column
    59         do jcol = istartcol,iendcol
     59        DO jcol = istartcol,iendcol
    6060
    6161          ! Compute relative humidity with respect to liquid
     
    7171          od_aerosol_mono = 0.0_jprb
    7272
    73           do jtype = 1,config%n_aerosol_types
     73          DO jtype = 1,config%n_aerosol_types
    7474            ! Add the optical depth for this aerosol type to the total for all aerosols. 
    7575            ! Note that the following expressions are array-wise
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_lw.F90

    r4773 r5158  
    8383    ! also the "source", which is the upwelling flux due to emission
    8484    ! below that level
    85     do jlev = nlev,1,-1
     85    DO jlev = nlev,1,-1
    8686      ! Next loop over columns. We could do this by indexing the
    8787      ! entire inner dimension as follows, e.g. for the first line:
     
    9090      ! routine by a factor of 2!  Rather, we do it with an explicit
    9191      ! loop.
    92       do jcol = 1,ncol
     92      DO jcol = 1,ncol
    9393        ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10:
    9494        inv_denominator(jcol,jlev) = 1.0_jprb &
     
    114114    ! Work back down through the atmosphere computing the fluxes at
    115115    ! each half-level
    116     do jlev = 1,nlev
    117       do jcol = 1,ncol
     116    DO jlev = 1,nlev
     117      DO jcol = 1,ncol
    118118        ! Shonk & Hogan (2008) Eq 14 (after simplification):
    119119        flux_dn(jcol,jlev+1) &
     
    201201    ! also the "source", which is the upwelling flux due to emission
    202202    ! below that level
    203     do jlev = nlev,i_cloud_top,-1
     203    DO jlev = nlev,i_cloud_top,-1
    204204      if (is_clear_sky_layer(jlev)) then
    205205        ! Reflectance of this layer is zero, simplifying the expression
    206         do jcol = 1,ncol
     206        DO jcol = 1,ncol
    207207          albedo(jcol,jlev) = transmittance(jcol,jlev)*transmittance(jcol,jlev)*albedo(jcol,jlev+1)
    208208          source(jcol,jlev) = source_up(jcol,jlev) &
     
    212212      else
    213213        ! Loop over columns; explicit loop seems to be faster
    214         do jcol = 1,ncol
     214        DO jcol = 1,ncol
    215215          ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10:
    216216          inv_denominator(jcol,jlev) = 1.0_jprb &
     
    231231    flux_up(:,i_cloud_top) = source(:,i_cloud_top) &
    232232         &                 + albedo(:,i_cloud_top)*flux_dn(:,i_cloud_top)
    233     do jlev = i_cloud_top-1,1,-1
     233    DO jlev = i_cloud_top-1,1,-1
    234234      flux_up(:,jlev) = transmittance(:,jlev)*flux_up(:,jlev+1) + source_up(:,jlev)
    235235    end do
     
    237237    ! Work back down through the atmosphere from cloud top computing
    238238    ! the fluxes at each half-level
    239     do jlev = i_cloud_top,nlev
     239    DO jlev = i_cloud_top,nlev
    240240      if (is_clear_sky_layer(jlev)) then
    241         do jcol = 1,ncol
     241        DO jcol = 1,ncol
    242242          flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) &
    243243               &               + source_dn(jcol,jlev)
     
    246246        end do
    247247      else
    248         do jcol = 1,ncol
     248        DO jcol = 1,ncol
    249249          ! Shonk & Hogan (2008) Eq 14 (after simplification):
    250250          flux_dn(jcol,jlev+1) &
     
    309309! Added for DWD (2020)
    310310!NEC$ outerloop_unroll(8)
    311     do jlev = 1,nlev
    312       do jcol = 1,ncol
     311    DO jlev = 1,nlev
     312      DO jcol = 1,ncol
    313313        flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) + source_dn(jcol,jlev)
    314314      end do
     
    322322! Added for DWD (2020)
    323323!NEC$ outerloop_unroll(8)
    324     do jlev = nlev,1,-1
    325       do jcol = 1,ncol
     324    DO jlev = nlev,1,-1
     325      DO jcol = 1,ncol
    326326        flux_up(jcol,jlev) = transmittance(jcol,jlev)*flux_up(jcol,jlev+1) + source_up(jcol,jlev)
    327327      end do
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_sw.F90

    r4773 r5158  
    8484    ! half-level by working down through the atmosphere
    8585    flux_dn_direct(:,1) = incoming_toa
    86     do jlev = 1,nlev
     86    DO jlev = 1,nlev
    8787      flux_dn_direct(:,jlev+1) = flux_dn_direct(:,jlev)*trans_dir_dir(:,jlev)
    8888    end do
     
    100100! Added for DWD (2020)
    101101!NEC$ outerloop_unroll(8)
    102     do jlev = nlev,1,-1
     102    DO jlev = nlev,1,-1
    103103      ! Next loop over columns. We could do this by indexing the
    104104      ! entire inner dimension as follows, e.g. for the first line:
     
    107107      ! routine by a factor of 2!  Rather, we do it with an explicit
    108108      ! loop.
    109       do jcol = 1,ncol
     109      DO jcol = 1,ncol
    110110        ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10:
    111111        inv_denominator(jcol,jlev) = 1.0_jprb / (1.0_jprb-albedo(jcol,jlev+1)*reflectance(jcol,jlev))
     
    132132! Added for DWD (2020)
    133133!NEC$ outerloop_unroll(8)
    134     do jlev = 1,nlev
    135       do jcol = 1,ncol
     134    DO jlev = 1,nlev
     135      DO jcol = 1,ncol
    136136        ! Shonk & Hogan (2008) Eq 14 (after simplification):
    137137        flux_dn_diffuse(jcol,jlev+1) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics.F90

    r4853 r5158  
    243243
    244244        if (ao%use_hydrophilic) then
    245           do jtype = 1,n_type_philic
     245          DO jtype = 1,n_type_philic
    246246            ao%mass_ext_sw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype))
    247247            ao%ssa_sw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype) &
     
    266266
    267267        if (ao%use_hydrophilic) then
    268           do jtype = 1,n_type_philic
     268          DO jtype = 1,n_type_philic
    269269            ao%mass_ext_lw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype))
    270270            ao%ssa_lw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype) &
     
    280280      if (allocated(ao%wavelength_mono)) then
    281281        ! Monochromatic wavelengths also required
    282         do jwl = 1,nmono
     282        DO jwl = 1,nmono
    283283          ! Wavelength (m) to wavenumber (cm-1)
    284284          wavenumber_target = 0.01_jprb / ao%wavelength_mono(jwl)
     
    292292          else
    293293            iwn = 1
    294             do while (wavenumber(iwn+1) < wavenumber_target .and. iwn < nwn-1)
     294            DO while (wavenumber(iwn+1) < wavenumber_target .and. iwn < nwn-1)
    295295              iwn = iwn + 1
    296296            end do
     
    419419
    420420      if (ao%use_hydrophilic) then
    421         do jtype = 1,ao%n_type_philic
     421        DO jtype = 1,ao%n_type_philic
    422422          ao%mass_ext_sw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_sw_philic(:,:,jtype))
    423423          ao%ssa_sw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_sw_philic(:,:,jtype) &
     
    451451
    452452      if (ao%use_hydrophilic) then
    453         do jtype = 1,ao%n_type_philic
     453        DO jtype = 1,ao%n_type_philic
    454454          ao%mass_ext_lw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_lw_philic(:,:,jtype))
    455455          ao%ssa_lw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_lw_philic(:,:,jtype) &
     
    579579      ! Aerosol mixing ratios have been provided
    580580
    581       do jtype = 1,config%n_aerosol_types
     581      DO jtype = 1,config%n_aerosol_types
    582582        if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then
    583583          write(nulerr,'(a)') '*** Error: not all aerosol types are defined'
     
    612612
    613613      ! Loop over column
    614       do jcol = istartcol,iendcol
     614      DO jcol = istartcol,iendcol
    615615
    616616        ! Reset temporary arrays
     
    622622        scat_g_lw_aerosol = 0.0_jprb
    623623
    624         do jlev = istartlev,iendlev
     624        DO jlev = istartlev,iendlev
    625625          ! Compute relative humidity with respect to liquid
    626626          ! saturation and the index to the relative-humidity index of
     
    634634        end do
    635635
    636         do jtype = 1,config%n_aerosol_types
     636        DO jtype = 1,config%n_aerosol_types
    637637          itype = ao%itype(jtype)
    638638
     
    643643          ! dimension being spectral band.
    644644          if (ao%iclass(jtype) == IAerosolClassHydrophobic) then
    645             do jlev = istartlev,iendlev
     645            DO jlev = istartlev,iendlev
    646646              mixing_ratio = aerosol%mixing_ratio(jcol,jlev,jtype)
    647               do jband = 1,config%n_bands_sw
     647              DO jband = 1,config%n_bands_sw
    648648                local_od_sw = factor(jlev) * mixing_ratio &
    649649                     &  * ao%mass_ext_sw_phobic(jband,itype)
     
    656656              end do
    657657              if (config%do_lw_aerosol_scattering) then
    658                 do jband = 1,config%n_bands_lw
     658                DO jband = 1,config%n_bands_lw
    659659                  local_od_lw = factor(jlev) * mixing_ratio &
    660660                       &  * ao%mass_ext_lw_phobic(jband,itype)
     
    670670                ! weight the optical depth by the single scattering
    671671                ! co-albedo
    672                 do jband = 1,config%n_bands_lw
     672                DO jband = 1,config%n_bands_lw
    673673                  od_lw_aerosol(jband,jlev) = od_lw_aerosol(jband,jlev) &
    674674                       &  + factor(jlev) * mixing_ratio &
     
    682682            ! Hydrophilic aerosols require the look-up tables to
    683683            ! be indexed with irh
    684             do jlev = istartlev,iendlev
     684            DO jlev = istartlev,iendlev
    685685              mixing_ratio = aerosol%mixing_ratio(jcol,jlev,jtype)
    686686              irh = irhs(jlev)
    687               do jband = 1,config%n_bands_sw
     687              DO jband = 1,config%n_bands_sw
    688688                local_od_sw = factor(jlev) * mixing_ratio &
    689689                     &  * ao%mass_ext_sw_philic(jband,irh,itype)
     
    696696              end do
    697697              if (config%do_lw_aerosol_scattering) then
    698                 do jband = 1,config%n_bands_lw
     698                DO jband = 1,config%n_bands_lw
    699699                  local_od_lw = factor(jlev) * mixing_ratio &
    700700                       &  * ao%mass_ext_lw_philic(jband,irh,itype)
     
    710710                ! weight the optical depth by the single scattering
    711711                ! co-albedo
    712                 do jband = 1,config%n_bands_lw
     712                DO jband = 1,config%n_bands_lw
    713713                  od_lw_aerosol(jband,jlev) = od_lw_aerosol(jband,jlev) &
    714714                       &  + factor(jlev) * mixing_ratio &
     
    740740
    741741          ! We can assume the band and g-point indices are the same
    742           do jlev = 1,nlev
    743             do jg = 1,config%n_g_sw
     742          DO jlev = 1,nlev
     743            DO jg = 1,config%n_g_sw
    744744              local_scat = ssa_sw(jg,jlev,jcol)*od_sw(jg,jlev,jcol) + scat_sw_aerosol(jg,jlev)
    745745              od_sw(jg,jlev,jcol) = od_sw(jg,jlev,jcol) + od_sw_aerosol(jg,jlev)
     
    751751        else
    752752
    753           do jlev = 1,nlev
    754             do jg = 1,config%n_g_sw
     753          DO jlev = 1,nlev
     754            DO jg = 1,config%n_g_sw
    755755              ! Need to map between bands and g-points
    756756              iband = config%i_band_from_reordered_g_sw(jg)
     
    781781               &                             scat_lw_aerosol, scat_g_lw_aerosol)
    782782
    783           do jlev = istartlev,iendlev
    784             do jg = 1,config%n_g_lw
     783          DO jlev = istartlev,iendlev
     784            DO jg = 1,config%n_g_lw
    785785              iband = config%i_band_from_reordered_g_lw(jg)
    786786              local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband,jlev)
     
    802802          if (config%do_cloud_aerosol_per_lw_g_point) then
    803803            ! We can assume band and g-point indices are the same
    804             do jlev = istartlev,iendlev
    805               do jg = 1,config%n_g_lw
     804            DO jlev = istartlev,iendlev
     805              DO jg = 1,config%n_g_lw
    806806                od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) + od_lw_aerosol(jg,jlev)
    807807              end do
    808808            end do
    809809          else
    810             do jlev = istartlev,iendlev
    811               do jg = 1,config%n_g_lw
     810            DO jlev = istartlev,iendlev
     811              DO jg = 1,config%n_g_lw
    812812                od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) &
    813813                     &  + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg),jlev)
     
    901901
    902902      ! Loop over position
    903       do jcol = istartcol,iendcol
     903      DO jcol = istartcol,iendcol
    904904! Added for DWD (2020)
    905905!NEC$ forced_collapse
    906         do jlev = istartlev,iendlev
    907           do jb = 1,config%n_bands_sw
     906        DO jlev = istartlev,iendlev
     907          DO jb = 1,config%n_bands_sw
    908908            od_sw_aerosol(jb,jlev) = aerosol%od_sw(jb,jlev,jcol)
    909909            scat_sw_aerosol(jb,jlev) = aerosol%ssa_sw(jb,jlev,jcol) * od_sw_aerosol(jb,jlev)
     
    923923        ! properties (noting that any gas scattering will have an
    924924        ! asymmetry factor of zero)
    925         do jlev = istartlev,iendlev
     925        DO jlev = istartlev,iendlev
    926926          if (od_sw_aerosol(1,jlev) > 0.0_jprb) then
    927             do jg = 1,config%n_g_sw
     927            DO jg = 1,config%n_g_sw
    928928              iband = config%i_band_from_reordered_g_sw(jg)
    929929              local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband,jlev)
     
    961961
    962962        ! Loop over position
    963         do jcol = istartcol,iendcol
     963        DO jcol = istartcol,iendcol
    964964! Added for DWD (2020)
    965965!NEC$ forced_collapse
    966           do jlev = istartlev,iendlev
    967             do jb = 1,config%n_bands_lw
     966          DO jlev = istartlev,iendlev
     967            DO jb = 1,config%n_bands_lw
    968968              od_lw_aerosol(jb,jlev) = aerosol%od_lw(jb,jlev,jcol)
    969969              scat_lw_aerosol(jb,jlev) = aerosol%ssa_lw(jb,jlev,jcol) * od_lw_aerosol(jb,jlev)
     
    974974            end do
    975975          end do
    976           do jlev = istartlev,iendlev
    977             do jg = 1,config%n_g_lw
     976          DO jlev = istartlev,iendlev
     977            DO jg = 1,config%n_g_lw
    978978              iband = config%i_band_from_reordered_g_lw(jg)
    979979              if (od_lw_aerosol(iband,jlev) > 0.0_jprb) then
     
    995995
    996996        ! Loop over position
    997         do jcol = istartcol,iendcol
     997        DO jcol = istartcol,iendcol
    998998! Added for DWD (2020)
    999999!NEC$ forced_collapse
    1000           do jlev = istartlev,iendlev
     1000          DO jlev = istartlev,iendlev
    10011001            ! If aerosol longwave scattering is not included then we
    10021002            ! weight the optical depth by the single scattering
    10031003            ! co-albedo
    1004             do jb = 1, config%n_bands_lw
     1004            DO jb = 1, config%n_bands_lw
    10051005              od_lw_aerosol(jb,jlev) = aerosol%od_lw(jb,jlev,jcol) &
    10061006                 &  * (1.0_jprb - aerosol%ssa_lw(jb,jlev,jcol))
    10071007            end do
    10081008          end do
    1009           do jlev = istartlev,iendlev
    1010             do jg = 1,config%n_g_lw
     1009          DO jlev = istartlev,iendlev
     1010            DO jg = 1,config%n_g_lw
    10111011              od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) &
    10121012                   &  + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg),jlev)
     
    11201120    if (lhook) call dr_hook('radiation_aerosol_optics:aerosol_extinction',0,hook_handle)
    11211121
    1122     do jtype = 1,config%n_aerosol_types
     1122    DO jtype = 1,config%n_aerosol_types
    11231123      if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then
    11241124        write(nulerr,'(a)') '*** Error: not all aerosol types are defined'
     
    11381138
    11391139    ! Loop over position
    1140     do jcol = istartcol,iendcol
     1140    DO jcol = istartcol,iendcol
    11411141      ext = 0.0_jprb
    11421142      ! Get relative-humidity index
    11431143      irh = ao%calc_rh_index(relative_humidity(jcol))
    11441144      ! Add extinction coefficients from each aerosol type
    1145       do jtype = 1,config%n_aerosol_types
     1145      DO jtype = 1,config%n_aerosol_types
    11461146        if (ao%iclass(jtype) == IAerosolClassHydrophobic) then
    11471147          ext = ext + mixing_ratio(jcol,jtype) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_data.F90

    r4860 r5158  
    619619    iend   = ubound(itypes,1)
    620620
    621     do jtype = istart, iend
     621    DO jtype = istart, iend
    622622      if (itypes(jtype) == 0) then
    623623        call this%set_empty_type(jtype)
     
    655655    else
    656656      calc_rh_index = 1
    657       do while (rh > this%rh_lower(calc_rh_index + 1))
     657      DO while (rh > this%rh_lower(calc_rh_index + 1))
    658658        calc_rh_index = calc_rh_index + 1
    659659      end do
     
    682682    end if
    683683
    684     do jtype = 1,size(i_type_map)
     684    DO jtype = 1,size(i_type_map)
    685685      if (i_type_map(jtype) > 0) then
    686686        write(nulout,'(i4,a,a)') jtype, ' -> hydrophobic type ', &
     
    714714
    715715    ! Find index of first character
    716     do while (i_line_current < iline)
     716    DO while (i_line_current < iline)
    717717      i_start_new = scan(str(istart:iend), new_line(' '))
    718718      if (i_start_new == 0) then
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_description.F90

    r4853 r5158  
    165165   
    166166    ! Loop over hydrophilic types
    167     do ja = 1,size(this%bin_philic)
     167    DO ja = 1,size(this%bin_philic)
    168168      ! Check if we have a match
    169169      if (to_string(this%code_philic(:,ja)) == code_str &
     
    176176    end do
    177177    ! Repeat for the hydrophobic types
    178     do ja = 1,size(this%bin_phobic)
     178    DO ja = 1,size(this%bin_phobic)
    179179      if (to_string(this%code_phobic(:,ja)) == code_str &
    180180           &  .and. trim(to_string(this%optical_model_phobic(:,ja))) &
     
    255255    if (lhydrophilic) then
    256256      ! Loop over hydrophilic aerosol types
    257       do ja = 1,size(this%bin_philic)
     257      DO ja = 1,size(this%bin_philic)
    258258        current_score = 0
    259259        if (to_string(this%code_philic(:,ja)) == code_str) then
     
    306306    else
    307307      ! Loop over hydrophobic aerosol types
    308       do ja = 1,size(this%bin_phobic)
     308      DO ja = 1,size(this%bin_phobic)
    309309        current_score = 0
    310310        if (to_string(this%code_phobic(:,ja)) == code_str) then
     
    376376    character(len=size(arr)) :: str
    377377    integer :: jc
    378     do jc = 1,size(arr)
     378    DO jc = 1,size(arr)
    379379      if (ichar(arr(jc)) == 0) then
    380380        ! Replace NULL character with a space
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud.F90

    r4911 r5158  
    251251      ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
    252252      ! don't take the logarithm of the first pressure in each column.
    253       do jcol = i1,i2
     253      DO jcol = i1,i2
    254254        this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length) &
    255255             &                            * thermodynamics%temperature_hl(jcol,2) &
     
    258258      end do
    259259
    260       do jlev = 2,nlev-1
    261         do jcol = i1,i2
     260      DO jlev = 2,nlev-1
     261        DO jcol = i1,i2
    262262          this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) &
    263263              &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
     
    271271       ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
    272272       ! don't take the logarithm of the last pressure in each column.
    273       do jlev = 1,nlev-2
    274         do jcol = i1,i2
     273      DO jlev = 1,nlev-2
     274        DO jcol = i1,i2
    275275          this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) &
    276276              &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
     
    280280      end do
    281281
    282       do jcol = i1,i2
     282      DO jcol = i1,i2
    283283        this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length) &
    284284            &                            * thermodynamics%temperature_hl(jcol,nlev) &
     
    340340      ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
    341341      ! don't take the logarithm of the first pressure in each column.
    342       do jcol = istartcol,iendcol
     342      DO jcol = istartcol,iendcol
    343343        this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length(jcol)) &
    344344             &                            * thermodynamics%temperature_hl(jcol,2) &
     
    347347      end do
    348348
    349       do jlev = 2,nlev-1
    350         do jcol = istartcol,iendcol
     349      DO jlev = 2,nlev-1
     350        DO jcol = istartcol,iendcol
    351351          this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol)) &
    352352              &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
     
    360360       ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
    361361       ! don't take the logarithm of the last pressure in each column.
    362       do jlev = 1,nlev-2
    363         do jcol = istartcol,iendcol
     362      DO jlev = 1,nlev-2
     363        DO jcol = istartcol,iendcol
    364364          this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol)) &
    365365              &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
     
    369369      end do
    370370
    371       do jcol = istartcol,iendcol
     371      DO jcol = istartcol,iendcol
    372372        this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length(jcol)) &
    373373            &                            * thermodynamics%temperature_hl(jcol,nlev) &
     
    423423      ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we
    424424      ! don't take the logarithm of the first pressure in each column.
    425       do jcol = istartcol,iendcol
     425      DO jcol = istartcol,iendcol
    426426        this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length(jcol,1)) &
    427427             &                            * thermodynamics%temperature_hl(jcol,2) &
     
    430430      end do
    431431
    432       do jlev = 2,nlev-1
    433         do jcol = istartcol,iendcol
     432      DO jlev = 2,nlev-1
     433        DO jcol = istartcol,iendcol
    434434          this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol,jlev)) &
    435435              &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
     
    443443       ! to top-of-atmosphere).  In case pressure_hl(:,nlev+1)=0, we
    444444       ! don't take the logarithm of the last pressure in each column.
    445       do jlev = 1,nlev-2
    446         do jcol = istartcol,iendcol
     445      DO jlev = 1,nlev-2
     446        DO jcol = istartcol,iendcol
    447447          this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol,jlev)) &
    448448              &                            * thermodynamics%temperature_hl(jcol,jlev+1) &
     
    452452      end do
    453453
    454       do jcol = istartcol,iendcol
     454      DO jcol = istartcol,iendcol
    455455      ! AI ATTENTION a verifier decorrelation_length(jcol,nlev-1) ou nlev
    456456        this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length(jcol,nlev-1)) &
     
    663663    end if
    664664
    665     do jcol = i1,i2
     665    DO jcol = i1,i2
    666666      eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) &
    667667           &  * (0.5_jprb / pressure_hl(jcol,isurf))
     
    761761    end if
    762762
    763     do jcol = i1,i2
     763    DO jcol = i1,i2
    764764      eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) &
    765765           &  * (0.5_jprb / pressure_hl(jcol,isurf))
     
    804804    ntype = size(this%mixing_ratio,3)
    805805   
    806     do jlev = 1,nlev
    807       do jcol = istartcol,iendcol
     806    DO jlev = 1,nlev
     807      DO jcol = istartcol,iendcol
    808808        sum_mixing_ratio(jcol) = 0.0_jprb
    809809      end do
    810       do jh = 1, ntype
    811         do jcol = istartcol,iendcol
     810      DO jh = 1, ntype
     811        DO jcol = istartcol,iendcol
    812812          sum_mixing_ratio(jcol) = sum_mixing_ratio(jcol) + this%mixing_ratio(jcol,jlev,jh)
    813813        end do
    814814      end do
    815       do jcol = istartcol,iendcol
     815      DO jcol = istartcol,iendcol
    816816        if (this%fraction(jcol,jlev)        < cloud_fraction_threshold &
    817817             &  .or. sum_mixing_ratio(jcol) < cloud_mixing_ratio_threshold) then
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_cover.F90

    r4853 r5158  
    159159    cum_product = 1.0_jprb - frac(1)
    160160    cum_cloud_cover(1) = frac(1)
    161     do jlev = 1,nlev-1
     161    DO jlev = 1,nlev-1
    162162      ! Compute the combined cloud cover of layers jlev and jlev+1
    163163      pair_cloud_cover(jlev) = max(frac(jlev),frac(jlev+1))
     
    242242    cum_product = 1.0_jprb - frac(1)
    243243    cum_cloud_cover(1) = frac(1)
    244     do jlev = 1,nlev-1
     244    DO jlev = 1,nlev-1
    245245      ! Convert to "alpha" overlap parameter if necessary
    246246      if (do_overlap_conversion) then
     
    258258#ifdef DWD_VECTOR_OPTIMIZATIONS
    259259    end do
    260     do jlev = 1,nlev-1
     260    DO jlev = 1,nlev-1
    261261#endif
    262262      if (frac(jlev) >= MaxCloudFrac) then
     
    385385    jlev = 1
    386386    nobj = 0
    387     do while (jlev <= nlev)
     387    DO while (jlev <= nlev)
    388388      if (frac(jlev) > min_frac) then
    389389        ! Starting a new object: store its top
     
    392392        ! Find its maximum cloud fraction
    393393        jlev = jlev + 1
    394         do while (jlev <= nlev)
     394        DO while (jlev <= nlev)
    395395          if (frac(jlev) < frac(jlev-1)) then
    396396            exit
     
    400400        i_max_obj(nobj) = jlev - 1
    401401        ! Find its base
    402         do while (jlev <= nlev)
     402        DO while (jlev <= nlev)
    403403          if (frac(jlev) > frac(jlev-1) .or. frac(jlev) <= min_frac) then
    404404            exit
     
    435435
    436436        ! Compute the combined cloud cover of pairs of layers
    437         do jlev = 1,nlev-1
     437        DO jlev = 1,nlev-1
    438438          pair_cloud_cover(jlev) &
    439439               &  = overlap_param(jlev)*max(frac(jlev),frac(jlev+1)) &
     
    444444        ! adjacent objects as the product of the layerwise overlap
    445445        ! parameters between their layers of maximum cloud fraction
    446         do jobj = 1,nobj-1
     446        DO jobj = 1,nobj-1
    447447          alpha_obj(jobj) &
    448448               &  = product(overlap_param(i_max_obj(jobj):i_max_obj(jobj+1)-1))
     
    456456             &                     frac(1:nlev-1), frac(2:nlev))
    457457        ! Compute the combined cloud cover of pairs of layers
    458         do jlev = 1,nlev-1
     458        DO jlev = 1,nlev-1
    459459          pair_cloud_cover(jlev) &
    460460               &  = overlap_alpha(jlev)*max(frac(jlev),frac(jlev+1)) &
     
    465465        ! adjacent objects as the product of the layerwise overlap
    466466        ! parameters between their layers of maximum cloud fraction
    467         do jobj = 1,nobj-1
     467        DO jobj = 1,nobj-1
    468468          alpha_obj(jobj) &
    469469               &  = product(overlap_alpha(i_max_obj(jobj):i_max_obj(jobj+1)-1))
     
    476476      ! of each object: this will later be converted to the cumulative
    477477      ! cloud cover working down from TOA
    478       do jobj = 1,nobj
     478      DO jobj = 1,nobj
    479479        cum_cloud_cover(i_top_obj(jobj)) = frac(i_top_obj(jobj))
    480         do jlev = i_top_obj(jobj), i_base_obj(jobj)-1
     480        DO jlev = i_top_obj(jobj), i_base_obj(jobj)-1
    481481          if (frac(jlev) >= MaxCloudFrac) then
    482482            ! Cloud cover has reached one
     
    496496      ! Sequentially combine objects until there is only one left
    497497      ! covering the full vertical extent of clouds in the profile
    498       do while (nobj > 1)
     498      DO while (nobj > 1)
    499499        ! Find the most correlated adjacent pair of objects
    500500        alpha_max = 0.0_jprb
     
    506506
    507507        jobj = 1
    508         do while (jobj < nobj)
     508        DO while (jobj < nobj)
    509509          if (alpha_obj(jobj) > alpha_max) then
    510510            alpha_max = alpha_obj(jobj)
     
    533533        ! Scale the combined cloud cover of the lower object to
    534534        ! account for its overlap with the upper object
    535         do jlev = i_top_obj(iobj2),i_base_obj(iobj2)
     535        DO jlev = i_top_obj(iobj2),i_base_obj(iobj2)
    536536          cum_cloud_cover(jlev) = cum_cloud_cover(i_base_obj(iobj1)) &
    537537               +  cum_cloud_cover(jlev) * scaling
     
    554554      ! Ensure that the combined cloud cover of pairs of layers is
    555555      ! consistent with the overhang
    556       do jlev = 1,nlev-1
     556      DO jlev = 1,nlev-1
    557557        pair_cloud_cover(jlev) = max(pair_cloud_cover(jlev), &
    558558             &     frac(jlev)+cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_generator.F90

    r4853 r5158  
    152152
    153153    total_cloud_cover = cum_cloud_cover(nlev);
    154     do jlev = 1,nlev-1
     154    DO jlev = 1,nlev-1
    155155      overhang(jlev) = cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev)
    156156    end do
     
    166166      ! Find range of cloudy layers
    167167      jlev = 1
    168       do while (frac(jlev) <= 0.0_jprb)
     168      DO while (frac(jlev) <= 0.0_jprb)
    169169        jlev = jlev + 1
    170170      end do
    171171      ibegin = jlev
    172172      iend = jlev
    173       do jlev = jlev+1,nlev
     173      DO jlev = jlev+1,nlev
    174174        if (frac(jlev) > 0.0_jprb) then
    175175          iend = jlev
     
    180180      overlap_param_inhom = overlap_param
    181181
    182       do jlev = ibegin,iend-1
     182      DO jlev = ibegin,iend-1
    183183        if (overlap_param(jlev) > 0.0_jprb) then
    184184          overlap_param_inhom(jlev) &
     
    208208       
    209209        ! Loop over ng columns
    210         do jg = 1,ng
     210        DO jg = 1,ng
    211211          ! Find the cloud top height corresponding to the current
    212212          ! random number, and store in itrigger
    213213          trigger = rand_top(jg) * total_cloud_cover
    214214          jlev = ibegin
    215           do while (trigger > cum_cloud_cover(jlev) .and. jlev < iend)
     215          DO while (trigger > cum_cloud_cover(jlev) .and. jlev < iend)
    216216            jlev = jlev + 1
    217217          end do
     
    328328    ! Loop from the layer below the local cloud top down to the
    329329    ! bottom-most cloudy layer
    330     do jlev = itrigger+1,iend+1
     330    DO jlev = itrigger+1,iend+1
    331331      do_fill_od_scaling = .false.
    332332      if (jlev <= iend) then
     
    367367
    368368        ! Loop through the sequence of cloudy layers
    369         do jcloud = 2,n_layers_to_scale
     369        DO jcloud = 2,n_layers_to_scale
    370370          ! Use second random number, and inhomogeneity overlap
    371371          ! parameter, to decide whether the first random number
     
    463463    ! Loop from the layer below the local cloud top down to the
    464464    ! bottom-most cloudy layer
    465     do jlev = itrigger+1,iend
     465    DO jlev = itrigger+1,iend
    466466      iy = iy+1
    467467      if (is_cloudy(jlev-1)) then
     
    494494       
    495495    ! Loop through the sequence of cloudy layers
    496     do jcloud = 2,n_layers_to_scale
     496    DO jcloud = 2,n_layers_to_scale
    497497      ! Use second random number, and inhomogeneity overlap
    498498      ! parameter, to decide whether the first random number
     
    679679
    680680    ! Loop down through layers starting at the first cloudy layer
    681     do jlev = ibegin,iend
     681    DO jlev = ibegin,iend
    682682
    683683      if (is_any_cloud(jlev)) then
     
    685685! Added for DWD (2020)
    686686!NEC$ shortloop
    687         do jg = 1,ng
     687        DO jg = 1,ng
    688688          ! The intention is that all these operations are vectorizable,
    689689          ! since all are vector operations on vectors of length ng...
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_optics.F90

    r4773 r5158  
    288288      end if
    289289
    290       do jlev = 1,nlev
    291         do jcol = istartcol,iendcol
     290      DO jlev = 1,nlev
     291        DO jcol = istartcol,iendcol
    292292          ! Only do anything if cloud is present (assume that
    293293          ! cloud%crop_cloud_fraction has already been called)
     
    459459  ! Added for DWD (2020)
    460460  !NEC$ shortloop
    461               do jb = 1, config%n_bands_lw
     461              DO jb = 1, config%n_bands_lw
    462462                od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) + od_lw_ice(jb)
    463463                if (scat_od_lw_liq(jb)+scat_od_lw_ice(jb) > 0.0_jprb) then
     
    477477  ! Added for DWD (2020)
    478478  !NEC$ shortloop
    479               do jb = 1, config%n_bands_lw
     479              DO jb = 1, config%n_bands_lw
    480480                od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) - scat_od_lw_liq(jb) &
    481481                      &                   + od_lw_ice(jb) - scat_od_lw_ice(jb)
     
    484484  ! Added for DWD (2020)
    485485  !NEC$ shortloop
    486             do jb = 1, config%n_bands_sw
     486            DO jb = 1, config%n_bands_sw
    487487              od_sw_cloud(jb,jlev,jcol) = od_sw_liq(jb) + od_sw_ice(jb)
    488488              g_sw_cloud(jb,jlev,jcol) = (g_sw_liq(jb) * scat_od_sw_liq(jb) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_lw.F90

    r4773 r5158  
    9595
    9696    ! Loop through columns
    97     do jcol = istartcol,iendcol
     97    DO jcol = istartcol,iendcol
    9898
    9999      ! Compute the reflectance and transmittance of all layers,
    100100      ! neglecting clouds
    101       do jlev = 1,nlev
     101      DO jlev = 1,nlev
    102102        if (config%do_lw_aerosol_scattering) then
    103103          ! Scattering case: first compute clear-sky reflectance,
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_sw.F90

    r4773 r5158  
    103103
    104104    ! Loop through columns
    105     do jcol = istartcol,iendcol
     105    DO jcol = istartcol,iendcol
    106106      ! Only perform calculation if sun above the horizon
    107107      if (single_level%cos_sza(jcol) > 0.0_jprb) then
     
    114114          ! Delta-Eddington scaling has already been performed to the
    115115          ! aerosol part of od, ssa and g
    116           do jlev = 1,nlev
     116          DO jlev = 1,nlev
    117117            call calc_two_stream_gammas_sw(ng, cos_sza, &
    118118                 &  ssa(:,jlev,jcol), g(:,jlev,jcol), &
     
    128128        else
    129129          ! Apply delta-Eddington scaling to the aerosol-gas mixture
    130           do jlev = 1,nlev
     130          DO jlev = 1,nlev
    131131            od_total  =  od(:,jlev,jcol)
    132132            ssa_total = ssa(:,jlev,jcol)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_config.F90

    r4853 r5158  
    16571657    ! contains the weights of interest.  We now find the non-zero weights
    16581658    nweights = 0
    1659     do jband = 1,size(mapping,2)
     1659    DO jband = 1,size(mapping,2)
    16601660      if (mapping(2,jband) > 0.0_jprb) then
    16611661        nweights = nweights+1
     
    16741674           &  wavenumber2, ' cm-1):'
    16751675      if (this%do_cloud_aerosol_per_sw_g_point) then
    1676         do jband = 1, nweights
     1676        DO jband = 1, nweights
    16771677          write(nulout, '(a,i0,a,f8.4)') '  Shortwave g point ', iband(jband), ': ', weight(jband)
    16781678        end do
    16791679      else
    1680         do jband = 1, nweights
     1680        DO jband = 1, nweights
    16811681          wavenumber1_band = this%gas_optics_sw%spectral_def%wavenumber1_band(iband(jband))
    16821682          wavenumber2_band = this%gas_optics_sw%spectral_def%wavenumber2_band(iband(jband))
     
    17281728    allocate(diag_ind(ninterval+2))
    17291729    diag_ind = 0
    1730     do jint = 1,ninterval+2
     1730    DO jint = 1,ninterval+2
    17311731      diag_ind(jint) = jint
    17321732    end do
     
    18941894    ! Count the number of albedo/emissivity intervals
    18951895    ninterval = 0
    1896     do jint = 1,NMaxAlbedoIntervals
     1896    DO jint = 1,NMaxAlbedoIntervals
    18971897      if (this%i_sw_albedo_index(jint) > 0) then
    18981898        ninterval = jint
     
    19481948        write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_albedo_from_band_sw), &
    19491949             &  ' shortwave intervals to albedo intervals:'
    1950         do jband = 1,size(this%i_albedo_from_band_sw)
     1950        DO jband = 1,size(this%i_albedo_from_band_sw)
    19511951          write(nulout,'(a,i0)',advance='no') ' ', this%i_albedo_from_band_sw(jband)
    19521952        end do
     
    19721972    ! Count the number of albedo/emissivity intervals
    19731973    ninterval = 0
    1974     do jint = 1,NMaxAlbedoIntervals
     1974    DO jint = 1,NMaxAlbedoIntervals
    19751975      if (this%i_lw_emiss_index(jint) > 0) then
    19761976        ninterval = jint
     
    20262026        write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_emiss_from_band_lw), &
    20272027             &  ' longwave intervals to emissivity intervals:'
    2028         do jband = 1,size(this%i_emiss_from_band_lw)
     2028        DO jband = 1,size(this%i_emiss_from_band_lw)
    20292029          write(nulout,'(a,i0)',advance='no') ' ', this%i_emiss_from_band_lw(jband)
    20302030        end do
     
    20562056      is_not_found = .true.
    20572057
    2058       do jc = 0,size(enum_str)-1
     2058      DO jc = 0,size(enum_str)-1
    20592059        if (trim(str) == trim(enum_str(jc))) then
    20602060          icode = jc
     
    20662066        write(nulerr,'(a,a,a,a,a)',advance='no') '*** Error: ', trim(var_name), &
    20672067             &  ' must be one of: "', enum_str(0), '"'
    2068         do jc = 1,size(enum_str)-1
     2068        DO jc = 1,size(enum_str)-1
    20692069          write(nulerr,'(a,a,a)',advance='no') ', "', trim(enum_str(jc)), '"'
    20702070        end do
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_delta_eddington.h

    r4773 r5158  
    8686  integer :: j
    8787
    88   do j = 1,ng
     88  DO j = 1,ng
    8989    g            = scat_od_g(j) / max(scat_od(j), 1.0e-24)
    9090    f            = g*g
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd.F90

    r4853 r5158  
    209209    istart = 1
    210210    this%i_gas_mapping = 0
    211     do jgas = 1, this%ngas
     211    DO jgas = 1, this%ngas
    212212      if (jgas < this%ngas) then
    213213        inext = istart + scan(constituent_id(istart:nchar), ' ')
     
    217217      ! Find gas code
    218218      i_gas_code = 0
    219       do jjgas = 1, NMaxGases
     219      DO jjgas = 1, NMaxGases
    220220        if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then
    221221          i_gas_code = jjgas
     
    264264         &  this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries'
    265265    write(nulout, '(a)') '  Gases:'
    266     do jgas = 1,this%ngas
     266    DO jgas = 1,this%ngas
    267267      if (this%single_gas(jgas)%i_gas_code > 0) then
    268268        write(nulout, '(a,a,a)', advance='no') '    ', &
     
    373373
    374374    ! Interpolate input SSI to regular wavenumber grid
    375     do jwav_grid = 1,nwav_grid
    376       do jwav = 1,nwav-1
     375    DO jwav_grid = 1,nwav_grid
     376      DO jwav = 1,nwav-1
    377377        if (wavenumber(jwav) < wavenumber_grid(jwav_grid) &
    378378             &  .and. wavenumber(jwav+1) >= wavenumber_grid(jwav_grid)) then
     
    426426           &  ReferenceTSI, ' W m-2, solar cycle amplitude (at solar maximum), update to original solar irradiance'
    427427      iband = 0
    428       do jg = 1,this%ng
     428      DO jg = 1,this%ng
    429429        if (this%spectral_def%i_band_number(jg) > iband) then
    430430          iband = this%spectral_def%i_band_number(jg)
     
    516516    global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass)
    517517
    518     do jcol = istartcol,iendcol
     518    DO jcol = istartcol,iendcol
    519519
    520520      log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)))
    521521
    522       do jlev = 1,nlev
     522      DO jlev = 1,nlev
    523523        ! Find interpolation points in pressure
    524524        pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) &
     
    546546      optical_depth_fl(:,:,jcol) = 0.0_jprb
    547547     
    548       do jgas = 1,this%ngas
     548      DO jgas = 1,this%ngas
    549549
    550550        associate (single_gas => this%single_gas(jgas))
     
    561561            end if
    562562           
    563             do jlev = 1,nlev
     563            DO jlev = 1,nlev
    564564              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
    565565                   &        + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
     
    581581            end if
    582582           
    583             do jlev = 1,nlev
     583            DO jlev = 1,nlev
    584584              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
    585585                   &        + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
     
    592592            ! Composite gases
    593593            molar_abs => this%single_gas(jgas)%molar_abs
    594             do jlev = 1,nlev
     594            DO jlev = 1,nlev
    595595              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
    596596                   &  + (simple_multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
     
    611611            molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
    612612            mole_frac1 = exp(single_gas%log_mole_frac1)
    613             do jlev = 1,nlev
     613            DO jlev = 1,nlev
    614614              ! Take care of mole_fraction == 0
    615615              log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode)*scaling, mole_frac1))
     
    628628              !      &                +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) &
    629629              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1)))
    630             do jlev = 1,nlev
     630            DO jlev = 1,nlev
    631631              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
    632632                   &  + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode) * scaling) * ( &
     
    651651      ! Rayleigh scattering
    652652      if (this%is_sw .and. present(rayleigh_od_fl)) then
    653         do jlev = 1,nlev
     653        DO jlev = 1,nlev
    654654          rayleigh_od_fl(:,jlev,jcol) = global_multiplier &
    655655               &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat
     
    729729    od_fl(:,:,:) = 0.0_jprb
    730730
    731     do jlev = 1,nlev
    732       do jcol = istartcol,iendcol
     731    DO jlev = 1,nlev
     732      DO jcol = istartcol,iendcol
    733733
    734734        log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,jlev)+pressure_hl(jcol,jlev+1)))
     
    758758    end do
    759759
    760     do jgas = 1,this%ngas
     760    DO jgas = 1,this%ngas
    761761
    762762      associate (single_gas => this%single_gas(jgas))
     
    768768          molar_abs => this%single_gas(jgas)%molar_abs
    769769
    770           do jlev = 1,nlev
    771             do jg = 1, this%ng
    772               do jcol = istartcol,iendcol
     770          DO jlev = 1,nlev
     771            DO jg = 1, this%ng
     772              DO jcol = istartcol,iendcol
    773773                multiplier = simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode)
    774774
     
    785785          molar_abs => this%single_gas(jgas)%molar_abs
    786786
    787           do jlev = 1,nlev
    788             do jg = 1, this%ng
    789               do jcol = istartcol,iendcol
     787          DO jlev = 1,nlev
     788            DO jg = 1, this%ng
     789              DO jcol = istartcol,iendcol
    790790                multiplier = simple_multiplier(jcol,jlev)  * (mole_fraction_fl(jcol,jlev,igascode) &
    791791                   &                            - single_gas%reference_mole_frac)
     
    804804          molar_abs => this%single_gas(jgas)%molar_abs
    805805
    806           do jlev = 1,nlev
    807             do jg = 1, this%ng
    808               do jcol = istartcol,iendcol
     806          DO jlev = 1,nlev
     807            DO jg = 1, this%ng
     808              DO jcol = istartcol,iendcol
    809809                od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
    810810                   &        + (simple_multiplier(jcol,jlev)*tw1(jcol,jlev)) *   &
     
    823823            mole_frac1 = exp(single_gas%log_mole_frac1)
    824824
    825             do jlev = 1,nlev
    826               do jcol = istartcol,iendcol
     825            DO jlev = 1,nlev
     826              DO jcol = istartcol,iendcol
    827827                ! Take care of mole_fraction == 0
    828828                log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1))
     
    835835            end do
    836836
    837           do jlev = 1,nlev
    838             do jg = 1, this%ng
     837          DO jlev = 1,nlev
     838            DO jg = 1, this%ng
    839839!NEC$ select_vector
    840               do jcol = istartcol,iendcol
     840              DO jcol = istartcol,iendcol
    841841
    842842                od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
     
    866866
    867867      ! Ensure the optical depth is not negative
    868       do jcol = istartcol,iendcol
    869         do jlev = 1,nlev
    870           do jg = 1, this%ng
     868      DO jcol = istartcol,iendcol
     869        DO jlev = 1,nlev
     870          DO jg = 1, this%ng
    871871            optical_depth_fl(jg,jlev,jcol) = max(0.0_jprb, od_fl(jcol,jg,jlev))
    872872          end do
     
    876876      ! Rayleigh scattering
    877877      if (this%is_sw .and. present(rayleigh_od_fl)) then
    878         do jcol = istartcol,iendcol
    879           do jlev = 1,nlev
    880             do jg = 1, this%ng
     878        DO jcol = istartcol,iendcol
     879          DO jlev = 1,nlev
     880            DO jg = 1, this%ng
    881881              rayleigh_od_fl(jg,jlev,jcol) = global_multiplier &
    882882               &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat(jg)
     
    908908    integer    :: it1, jt
    909909
    910     do jt = 1,nt
     910    DO jt = 1,nt
    911911      tindex1 = (temperature(jt) - this%temperature1_planck) &
    912912           &   * (1.0_jprb / this%d_temperature_planck)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd_interface.F90

    r4853 r5158  
    272272      ! rayleigh optical depth: convert to total optical depth and
    273273      ! single-scattering albedo
    274       do jcol = istartcol,iendcol
    275         do jlev = 1, nlev
    276           do jg = 1, config%n_g_sw
     274      DO jcol = istartcol,iendcol
     275        DO jlev = 1, nlev
     276          DO jg = 1, config%n_g_sw
    277277            od_sw(jg,jlev,jcol)  = od_sw(jg,jlev,jcol) + ssa_sw(jg,jlev,jcol)
    278278            ssa_sw(jg,jlev,jcol) = ssa_sw(jg,jlev,jcol) / od_sw(jg,jlev,jcol)
     
    308308
    309309      ! Calculate the Planck function for each g point
    310       do jcol = istartcol,iendcol
     310      DO jcol = istartcol,iendcol
    311311        call config%gas_optics_lw%calc_planck_function(nlev+1, &
    312312             &  thermodynamics%temperature_hl(jcol,:), planck_hl(:,:,jcol))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_flux.F90

    r4853 r5158  
    423423             &               config%i_band_from_reordered_g_sw, &
    424424             &               this%sw_dn_surf_band, istartcol, iendcol)
    425         do jcol = istartcol,iendcol
     425        DO jcol = istartcol,iendcol
    426426          this%sw_dn_surf_band(:,jcol) &
    427427               &  = this%sw_dn_surf_band(:,jcol) &
     
    429429        end do
    430430      else
    431         do jcol = istartcol,iendcol
     431        DO jcol = istartcol,iendcol
    432432          call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
    433433               &           config%i_band_from_reordered_g_sw, &
     
    450450               &               config%i_band_from_reordered_g_sw, &
    451451               &               this%sw_dn_surf_clear_band, istartcol, iendcol)
    452           do jcol = istartcol,iendcol
     452          DO jcol = istartcol,iendcol
    453453            this%sw_dn_surf_clear_band(:,jcol) &
    454454                 &  = this%sw_dn_surf_clear_band(:,jcol) &
     
    456456          end do
    457457        else
    458           do jcol = istartcol,iendcol
     458          DO jcol = istartcol,iendcol
    459459            call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), &
    460460                 &           config%i_band_from_reordered_g_sw, &
     
    486486               &               this%sw_dn_diffuse_surf_canopy, istartcol, iendcol)
    487487        else
    488           do jcol = istartcol,iendcol
     488          DO jcol = istartcol,iendcol
    489489            call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
    490490                 &           config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
     
    502502        this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
    503503        this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = 0.0_jprb
    504         do jband = 1,config%n_bands_sw
    505           do jalbedoband = 1,nalbedoband
     504        DO jband = 1,config%n_bands_sw
     505          DO jalbedoband = 1,nalbedoband
    506506            if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
    507507              ! Initially, "diffuse" is actually "total"
     
    534534               &               this%lw_dn_surf_canopy, istartcol, iendcol)
    535535        else
    536           do jcol = istartcol,iendcol
     536          DO jcol = istartcol,iendcol
    537537            call indexed_sum(this%lw_dn_surf_g(:,jcol), &
    538538                 &           config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
     
    548548               &               lw_dn_surf_band, istartcol, iendcol)
    549549        else
    550           do jcol = istartcol,iendcol
     550          DO jcol = istartcol,iendcol
    551551            call indexed_sum(this%lw_dn_surf_g(:,jcol), &
    552552                 &           config%i_band_from_reordered_g_lw, &
     
    556556        nalbedoband = size(config%lw_emiss_weights,1)
    557557        this%lw_dn_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
    558         do jband = 1,config%n_bands_lw
    559           do jalbedoband = 1,nalbedoband
     558        DO jband = 1,config%n_bands_lw
     559          DO jalbedoband = 1,nalbedoband
    560560            if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then
    561561              this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) &
     
    602602             &               this%sw_up_toa_band, istartcol, iendcol)
    603603      else
    604         do jcol = istartcol,iendcol
     604        DO jcol = istartcol,iendcol
    605605          call indexed_sum(this%sw_dn_toa_g(:,jcol), &
    606606               &           config%i_band_from_reordered_g_sw, &
     
    618618               &               this%sw_up_toa_clear_band, istartcol, iendcol)
    619619        else
    620           do jcol = istartcol,iendcol
     620          DO jcol = istartcol,iendcol
    621621            call indexed_sum(this%sw_up_toa_clear_g(:,jcol), &
    622622                 &               config%i_band_from_reordered_g_sw, &
     
    634634             &               this%lw_up_toa_band, istartcol, iendcol)
    635635      else
    636         do jcol = istartcol,iendcol
     636        DO jcol = istartcol,iendcol
    637637          call indexed_sum(this%lw_up_toa_g(:,jcol), &
    638638               &           config%i_band_from_reordered_g_lw, &
     
    647647               &               this%lw_up_toa_clear_band, istartcol, iendcol)
    648648        else
    649           do jcol = istartcol,iendcol
     649          DO jcol = istartcol,iendcol
    650650            call indexed_sum(this%lw_up_toa_clear_g(:,jcol), &
    651651                 &               config%i_band_from_reordered_g_lw, &
     
    753753    iend   = ubound(source,1)
    754754
    755     do jg = istart, iend
     755    DO jg = istart, iend
    756756      ig = ind(jg)
    757757      dest(ig) = dest(ig) + source(jg)
     
    777777    iend   = ubound(source,1)
    778778
    779     do jg = istart, iend
     779    DO jg = istart, iend
    780780      ig = ind(jg)
    781781      dest(ig) = dest(ig) + source(jg)
     
    797797    dest = 0.0
    798798
    799     do jg = lbound(source,1), ubound(source,1)
     799    DO jg = lbound(source,1), ubound(source,1)
    800800      ig = ind(jg)
    801       do jc = ist, iend
     801      DO jc = ist, iend
    802802        dest(ig,jc) = dest(ig,jc) + source(jg,jc)
    803803      end do
     
    820820    nlev   = size(source,2)
    821821
    822     do jlev = 1,nlev
    823       do jg = istart, iend
     822    DO jlev = 1,nlev
     823      DO jg = istart, iend
    824824        ig = ind(jg)
    825825        dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev)
     
    846846    nlev   = size(source,2)
    847847
    848     do jlev = 1,nlev
    849       do jg = istart, iend
     848    DO jlev = 1,nlev
     849      DO jg = istart, iend
    850850        ig = ind(jg)
    851851        dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_gas.F90

    r4853 r5158  
    209209    this%is_well_mixed(igas) = .false.
    210210
    211     do jk = 1,this%nlev
    212       do jc = i1,i2
     211    DO jk = 1,this%nlev
     212      DO jc = i1,i2
    213213        this%mixing_ratio(jc,jk,igas) = mixing_ratio(jc-i1+1,jk)
    214214      end do
     
    294294    this%is_well_mixed(igas)           = .true.
    295295
    296     do jk = 1,this%nlev
    297       do jc = i1,i2
     296    DO jk = 1,this%nlev
     297      DO jc = i1,i2
    298298        this%mixing_ratio(jc,jk,igas) = mixing_ratio
    299299      end do
     
    397397      end if
    398398    else
    399       do jg = 1,this%ntype
     399      DO jg = 1,this%ntype
    400400        call this%set_units(iunits, igas=this%icode(jg), scale_factor=new_sf)
    401401      end do
     
    416416   
    417417    scaling = this%scale_factor
    418     do jg = 1,NMaxGases
     418    DO jg = 1,NMaxGases
    419419      if (iunits == IMassMixingRatio .and. this%iunits(jg) == IVolumeMixingRatio) then
    420420        scaling(jg) = scaling(jg) * GasMolarMass(jg) / AirMolarMass
     
    481481      end if
    482482    else
    483       do jg = 1,this%ntype
     483      DO jg = 1,this%ntype
    484484        call this%assert_units(iunits, igas=this%icode(jg), scale_factor=sf, istatus=istatus)
    485485      end do
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics.F90

    r4853 r5158  
    5353    ! Count number of cloud types
    5454    config%n_cloud_types = 0
    55     do jtype = 1,NMaxCloudTypes
     55    DO jtype = 1,NMaxCloudTypes
    5656      if (len_trim(config%cloud_type_name(jtype)) > 0) then
    5757        config%n_cloud_types = jtype
     
    8383
    8484    ! Load cloud optics data
    85     do jtype = 1,config%n_cloud_types
     85    DO jtype = 1,config%n_cloud_types
    8686      if (config%cloud_type_name(jtype)(1:1) == '/') then
    8787        file_name = trim(config%cloud_type_name(jtype))
     
    191191
    192192    ! Loop over cloud types
    193     do jtype = 1,config%n_cloud_types
     193    DO jtype = 1,config%n_cloud_types
    194194      ! Compute in-cloud water path
    195195      if (config%is_homogeneous) then
     
    238238    ! Scale the combined longwave optical properties
    239239    if (config%do_lw_cloud_scattering) then
    240       do jcol = istartcol, iendcol
    241         do jlev = 1,nlev
     240      DO jcol = istartcol, iendcol
     241        DO jlev = 1,nlev
    242242          if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
    243243            ! Note that original cloud optics does not do
     
    259259    if (config%do_sw) then
    260260      if (.not. config%do_sw_delta_scaling_with_gases) then
    261         do jcol = istartcol, iendcol
    262           do jlev = 1,nlev
     261        DO jcol = istartcol, iendcol
     262          DO jlev = 1,nlev
    263263            if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
    264264              call delta_eddington_extensive(od_sw_cloud(:,jlev,jcol), &
     
    269269      end if
    270270
    271       do jcol = istartcol, iendcol
    272         do jlev = 1,nlev
     271      DO jcol = istartcol, iendcol
     272        DO jlev = 1,nlev
    273273          if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
    274274            ! Scale to get asymmetry factor and single scattering albedo
    275             do jg = 1, config%n_bands_sw
     275            DO jg = 1, config%n_bands_sw
    276276              g_sw_cloud(jg,jlev,jcol) = g_sw_cloud(jg,jlev,jcol) &
    277277                 &  / max(ssa_sw_cloud(jg,jlev,jcol), 1.0e-15_jprb)
     
    308308    if (lhook) call dr_hook('radiation_general_cloud_optics:save',0,hook_handle)
    309309
    310     do jtype = 1,config%n_cloud_types
     310    DO jtype = 1,config%n_cloud_types
    311311      if (config%do_sw) then
    312312        associate(co_sw => config%cloud_optics_sw(jtype))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics_data.F90

    r4936 r5158  
    285285
    286286    if (present(scat_od)) then
    287       do jcol = 1,ncol
    288         do jlev = 1,nlev
     287      DO jcol = 1,ncol
     288        DO jlev = 1,nlev
    289289          if (cloud_fraction(jcol, jlev) > 0.0_jprb) then
    290290            re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) &
     
    293293            weight2 = re_index - ire
    294294            weight1 = 1.0_jprb - weight2
    295             do jg = 1, ng
     295            DO jg = 1, ng
    296296              od_local = water_path(jcol, jlev) * (weight1*this%mass_ext(jg,ire) &
    297297                 &                                +weight2*this%mass_ext(jg,ire+1))
     
    310310    else
    311311      ! No scattering: return the absorption optical depth
    312       do jcol = 1,ncol
    313         do jlev = 1,nlev
     312      DO jcol = 1,ncol
     313        DO jlev = 1,nlev
    314314          if (water_path(jcol, jlev) > 0.0_jprb) then
    315315            re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) &
     
    404404
    405405    ! Define effective radius
    406     do ire = 1,this%n_effective_radius
     406    DO ire = 1,this%n_effective_radius
    407407      effective_radius(ire) = this%effective_radius_0 + this%d_effective_radius*(ire-1)
    408408    end do
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_lw.F90

    r4773 r5158  
    119119
    120120    ! Loop through columns
    121     do jcol = istartcol,iendcol
     121    DO jcol = istartcol,iendcol
    122122
    123123      ! Is there any cloud in the profile?
    124124      is_cloudy_profile = .false.
    125       do jlev = 1,nlev
     125      DO jlev = 1,nlev
    126126        if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
    127127          is_cloudy_profile = .true.
     
    135135      ! the clear-sky layers since these will be needed when we come
    136136      ! to do the total-sky fluxes.
    137       do jlev = 1,nlev
     137      DO jlev = 1,nlev
    138138        if (config%do_clear .or. cloud%fraction(jcol,jlev) &
    139139             &                 < config%cloud_fraction_threshold) then
     
    201201      ! now.
    202202      if (is_cloudy_profile .or. .not. config%do_clear) then
    203         do jlev = 1,nlev
     203        DO jlev = 1,nlev
    204204          ! Compute combined gas+aerosol+cloud optical properties;
    205205          ! note that for clear layers, the reflectance and
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_sw.F90

    r4773 r5158  
    123123
    124124    ! Loop through columns
    125     do jcol = istartcol,iendcol
     125    DO jcol = istartcol,iendcol
    126126      ! Only perform calculation if sun above the horizon
    127127      if (single_level%cos_sza(jcol) > 0.0_jprb) then
     
    131131        ! Is there any cloud in the profile?
    132132        is_cloudy_profile = .false.
    133         do jlev = 1,nlev
     133        DO jlev = 1,nlev
    134134          if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
    135135            is_cloudy_profile = .true.
     
    146146          ! Delta-Eddington scaling has already been performed to the
    147147          ! aerosol part of od, ssa and g
    148           do jlev = 1,nlev
     148          DO jlev = 1,nlev
    149149            if (config%do_clear .or. cloud%fraction(jcol,jlev) &
    150150                 &                 < config%cloud_fraction_threshold) then
     
    164164        else
    165165          ! Apply delta-Eddington scaling to the aerosol-gas mixture
    166           do jlev = 1,nlev
     166          DO jlev = 1,nlev
    167167            if (config%do_clear .or. cloud%fraction(jcol,jlev) &
    168168                 &                 < config%cloud_fraction_threshold) then
     
    229229        ! fluxes now.
    230230        if (is_cloudy_profile .or. .not. config%do_clear) then
    231           do jlev = 1,nlev
     231          DO jlev = 1,nlev
    232232            ! Compute combined gas+aerosol+cloud optical properties;
    233233            ! note that for clear layers, the reflectance and
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ice_optics_fu.F90

    r4773 r5158  
    7373! Added for DWD (2020)
    7474!NEC$ shortloop
    75     do jb = 1, nb
     75    DO jb = 1, nb
    7676      od(jb) = iwp_gm_2 * (coeff(jb,1) + coeff(jb,2) * inv_de_um)
    7777      scat_od(jb) = od(jb) * (1.0_jprb - (coeff(jb,3) + de_um*(coeff(jb,4) &
     
    124124! Added for DWD (2020)
    125125!NEC$ shortloop
    126     do jb = 1, nb
     126    DO jb = 1, nb
    127127      od(jb) = iwp_gm_2 * (coeff(jb,1) + inv_de_um*(coeff(jb,2) &
    128128         &  + inv_de_um*coeff(jb,3)))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ifs_rrtm.F90

    r4853 r5158  
    386386    ZONEMINUS_ARRAY = ZONEMINUS
    387387
    388     do jlev=1,nlev
    389       do jcol= istartcol,iendcol
     388    DO jlev=1,nlev
     389      DO jcol= istartcol,iendcol
    390390        pressure_fl(jcol,jlev) &
    391391            &  = 0.5_jprb * (thermodynamics%pressure_hl(jcol,jlev+istartlev-1) &
     
    486486        ! arrays have dimensions in a different order to the inputs,
    487487        ! so there is some inefficiency here.
    488         do jgreorder = 1,config%n_g_lw
     488        DO jgreorder = 1,config%n_g_lw
    489489          iband = config%i_band_from_reordered_g_lw(jgreorder)
    490490          ig = config%i_g_from_reordered_g_lw(jgreorder)
    491491         
    492492          ! Top-of-atmosphere half level
    493           do jlev = 1,nlev
    494             do jcol = istartcol,iendcol
     493          DO jlev = 1,nlev
     494            DO jcol = istartcol,iendcol
    495495              ! Some g points can return negative optical depths;
    496496              ! specifically original g points 54-56 which causes
     
    504504      else
    505505        ! G points have not been reordered
    506         do jcol = istartcol,iendcol
    507           do jlev = 1,nlev
     506        DO jcol = istartcol,iendcol
     507          DO jlev = 1,nlev
    508508            ! Check for negative optical depth
    509509            od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
     
    544544      ! Scale the incoming solar per band, if requested
    545545      if (config%use_spectral_solar_scaling) then
    546         do jg = 1,JPGPT_SW
    547           do jcol = istartcol,iendcol
     546        DO jg = 1,JPGPT_SW
     547          DO jcol = istartcol,iendcol
    548548            ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
    549549                 &   single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
     
    556556      ! ZINCSOL will be zero.
    557557      if (present(incoming_sw)) then
    558         do jcol = istartcol,iendcol
     558        DO jcol = istartcol,iendcol
    559559          if (single_level%cos_sza(jcol) > 0.0_jprb) then
    560560! Added for DWD (2020)
     
    570570!      if (.true.) then
    571571        ! Account for reordered g points
    572         do jgreorder = 1,config%n_g_sw
     572        DO jgreorder = 1,config%n_g_sw
    573573          ig = config%i_g_from_reordered_g_sw(jgreorder)
    574           do jlev = 1,nlev
    575             do jcol = istartcol,iendcol
     574          DO jlev = 1,nlev
     575            DO jcol = istartcol,iendcol
    576576              ! Check for negative optical depth
    577577              od_sw (jgreorder,nlev+1-jlev,jcol) &
     
    587587      else
    588588        ! G points have not been reordered
    589         do jcol = istartcol,iendcol
    590           do jlev = 1,nlev
    591             do jg = 1,config%n_g_sw
     589        DO jcol = istartcol,iendcol
     590          DO jlev = 1,nlev
     591            DO jg = 1,config%n_g_sw
    592592              ! Check for negative optical depth
    593593              od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
     
    596596          end do
    597597          if (present(incoming_sw)) then
    598             do jg = 1,config%n_g_sw
     598            DO jg = 1,config%n_g_sw
    599599              incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
    600600            end do
     
    667667    ! lowest interpolation bound, and the fraction into interpolation
    668668    ! interval
    669     do jlev = 1,nlev+1
    670       do jcol = istartcol,iendcol
     669    DO jlev = 1,nlev+1
     670      DO jcol = istartcol,iendcol
    671671        temperature = thermodynamics%temperature_hl(jcol,jlev+ilevoffset)
    672672        if (temperature < 339.0_jprb .and. temperature >= 160.0_jprb) then
     
    687687
    688688      ! Calculate Planck functions per band
    689       do jband = 1,config%n_bands_lw
     689      DO jband = 1,config%n_bands_lw
    690690        factor = zfluxfac * delwave(jband)
    691         do jcol = istartcol,iendcol
     691        DO jcol = istartcol,iendcol
    692692          planck_store(jcol,jband) = factor &
    693693               &  * (totplnk(ind(jcol),jband) &
     
    706706          ! Top-of-atmosphere half level - note that PFRAC is on model
    707707          ! levels not half levels
    708           do jgreorder = 1,config%n_g_lw
     708          DO jgreorder = 1,config%n_g_lw
    709709            iband = config%i_band_from_reordered_g_lw(jgreorder)
    710710            ig = config%i_g_from_reordered_g_lw(jgreorder)
     
    713713          end do
    714714        else
    715           do jgreorder = 1,config%n_g_lw
     715          DO jgreorder = 1,config%n_g_lw
    716716            iband = config%i_band_from_reordered_g_lw(jgreorder)
    717717            ig = config%i_g_from_reordered_g_lw(jgreorder)
     
    726726          ! Top-of-atmosphere half level - note that PFRAC is on model
    727727          ! levels not half levels
    728           do jg = 1,config%n_g_lw
     728          DO jg = 1,config%n_g_lw
    729729            iband = config%i_band_from_g_lw(jg)
    730730            planck_hl(jg,1,:) = planck_store(:,iband) * PFRAC(:,jg,nlev)
    731731          end do
    732732        else
    733           do jg = 1,config%n_g_lw
     733          DO jg = 1,config%n_g_lw
    734734            iband = config%i_band_from_g_lw(jg)
    735735            planck_tmp(:,jg) = planck_store(:,iband) * PFRAC(:,jg,nlev+2-jlev)
    736736          end do
    737           do jcol = istartcol,iendcol
     737          DO jcol = istartcol,iendcol
    738738            planck_hl(:,jlev,jcol) = planck_tmp(jcol,:)
    739739          end do
     
    794794
    795795    ! Work out surface interpolations
    796     do jcol = istartcol,iendcol
     796    DO jcol = istartcol,iendcol
    797797      Tsurf = temperature(jcol)
    798798      if (Tsurf < 339.0_jprb .and. Tsurf >= 160.0_jprb) then
     
    813813
    814814    ! Calculate Planck functions per band
    815     do jband = 1,config%n_bands_lw
     815    DO jband = 1,config%n_bands_lw
    816816      factor = zfluxfac * delwave(jband)
    817       do jcol = istartcol,iendcol
     817      DO jcol = istartcol,iendcol
    818818        planck_store(jcol,jband) = factor &
    819819             &  * (totplnk(ind(jcol),jband) &
     
    829829      ! in pressure since the the functions above treat pressure
    830830      ! decreasing with increasing index.
    831       do jgreorder = 1,config%n_g_lw
     831      DO jgreorder = 1,config%n_g_lw
    832832        iband = config%i_band_from_reordered_g_lw(jgreorder)
    833833        ig = config%i_g_from_reordered_g_lw(jgreorder)
     
    836836    else
    837837      ! G points have not been reordered
    838       do jg = 1,config%n_g_lw
     838      DO jg = 1,config%n_g_lw
    839839        iband = config%i_band_from_g_lw(jg)
    840840        planck_surf(jg,:) = planck_store(:,iband) * PFRAC(:,jg)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_liquid_optics_socrates.F90

    r4773 r5158  
    6565! Added for DWD (2020)
    6666!NEC$ shortloop
    67     do jb = 1, nb
     67    DO jb = 1, nb
    6868      od(jb) = lwp * (coeff(jb,1) + re*(coeff(jb,2) + re*coeff(jb,3))) &
    6969         &  / (1.0_jprb + re*(coeff(jb,4) + re*(coeff(jb,5) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_lw_derivatives.F90

    r4773 r5158  
    7474    ! Move up through the atmosphere computing the derivatives at each
    7575    ! half-level
    76     do jlev = nlev,1,-1
     76    DO jlev = nlev,1,-1
    7777      lw_derivatives_g = lw_derivatives_g * transmittance(:,jlev)
    7878      lw_derivatives(icol,jlev) = sum(lw_derivatives_g)
     
    122122    ! Move up through the atmosphere computing the derivatives at each
    123123    ! half-level
    124     do jlev = nlev,1,-1
     124    DO jlev = nlev,1,-1
    125125      lw_derivatives_g = lw_derivatives_g * transmittance(:,jlev)
    126126      lw_derivatives(icol,jlev) = (1.0_jprb - weight) * lw_derivatives(icol,jlev) &
     
    177177    ! Move up through the atmosphere computing the derivatives at each
    178178    ! half-level
    179     do jlev = nlev,1,-1
     179    DO jlev = nlev,1,-1
    180180      ! Compute effect of overlap at half-level jlev+1, yielding
    181181      ! derivatives just above that half-level
     
    243243      ! multiplication inline
    244244     
    245       do jlev = nlev,1,-1
     245      DO jlev = nlev,1,-1
    246246        ! Compute effect of overlap at half-level jlev+1, yielding
    247247        ! derivatives just above that half-level
     
    249249       
    250250        associate(A=>u_matrix(:,:,jlev+1), b=>lw_deriv_below)
    251           do jg = 1,ng   
     251          DO jg = 1,ng
    252252            ! Both inner and outer loop of the matrix loops j1 and j2 unrolled
    253253            ! inner loop:        j2=1             j2=2             j2=3
     
    273273      ! Move up through the atmosphere computing the derivatives at each
    274274      ! half-level
    275       do jlev = nlev,1,-1
     275      DO jlev = nlev,1,-1
    276276        ! Compute effect of overlap at half-level jlev+1, yielding
    277277        ! derivatives just above that half-level
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_matrix.F90

    r4773 r5158  
    9191      mat_x_vec(1:iend,1) = A(1:iend,1,1)*b(1:iend,1)
    9292    else
    93       do j1 = 1,m
    94         do j2 = 1,m
     93      DO j1 = 1,m
     94        DO j2 = 1,m
    9595          mat_x_vec(1:iend,j1) = mat_x_vec(1:iend,j1) &
    9696               &               + A(1:iend,j1,j2)*b(1:iend,j2)
     
    125125    singlemat_x_vec = 0.0_jprb
    126126
    127     do j1 = 1,m
    128       do j2 = 1,m
     127    DO j1 = 1,m
     128      DO j2 = 1,m
    129129        singlemat_x_vec(1:iend,j1) = singlemat_x_vec(1:iend,j1) &
    130130             &                    + A(j1,j2)*b(1:iend,j2)
     
    176176      m2block = 2*mblock
    177177      ! Do the top-left (C, D, F, G)
    178       do j2 = 1,m2block
    179         do j1 = 1,m2block
    180           do j3 = 1,m2block
     178      DO j2 = 1,m2block
     179        DO j1 = 1,m2block
     180          DO j3 = 1,m2block
    181181            mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) &
    182182                 &                  + A(1:iend,j1,j3)*B(1:iend,j3,j2)
     
    184184        end do
    185185      end do
    186       do j2 = m2block+1,m
     186      DO j2 = m2block+1,m
    187187        ! Do the top-right (E & H)
    188         do j1 = 1,m2block
    189           do j3 = 1,m
     188        DO j1 = 1,m2block
     189          DO j3 = 1,m
    190190            mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) &
    191191                 &                  + A(1:iend,j1,j3)*B(1:iend,j3,j2)
     
    193193        end do
    194194        ! Do the bottom-right (I)
    195         do j1 = m2block+1,m
    196           do j3 = m2block+1,m
     195        DO j1 = m2block+1,m
     196          DO j3 = m2block+1,m
    197197            mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) &
    198198                 &                  + A(1:iend,j1,j3)*B(1:iend,j3,j2)
     
    202202    else
    203203      ! Ordinary dense matrix
    204       do j2 = 1,m
    205         do j1 = 1,m
    206           do j3 = 1,m
     204      DO j2 = 1,m
     205        DO j1 = 1,m
     206          DO j3 = 1,m
    207207            mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) &
    208208                 &                  + A(1:iend,j1,j3)*B(1:iend,j3,j2)
     
    238238    singlemat_x_mat = 0.0_jprb
    239239
    240     do j2 = 1,m
    241       do j1 = 1,m
    242         do j3 = 1,m
     240    DO j2 = 1,m
     241      DO j1 = 1,m
     242        DO j3 = 1,m
    243243          singlemat_x_mat(1:iend,j1,j2) = singlemat_x_mat(1:iend,j1,j2) &
    244244               &                        + A(j1,j3)*B(1:iend,j3,j2)
     
    273273    mat_x_singlemat = 0.0_jprb
    274274
    275     do j2 = 1,m
    276       do j1 = 1,m
    277         do j3 = 1,m
     275    DO j2 = 1,m
     276      DO j1 = 1,m
     277        DO j3 = 1,m
    278278          mat_x_singlemat(1:iend,j1,j2) = mat_x_singlemat(1:iend,j1,j2) &
    279279               &                        + A(1:iend,j1,j3)*B(j3,j2)
     
    311311
    312312    identity_minus_mat_x_mat = - identity_minus_mat_x_mat
    313     do j = 1,m
     313    DO j = 1,m
    314314      identity_minus_mat_x_mat(1:iend,j,j) &
    315315           &     = 1.0_jprb + identity_minus_mat_x_mat(1:iend,j,j)
     
    337337   
    338338    sparse_x_dense = 0.0_jprb
    339     do j2 = 1,n2
    340       do j1 = 1,n1
     339    DO j2 = 1,n2
     340      DO j1 = 1,n1
    341341        if (sparse(j1,j2) /= 0.0_jprb) then
    342342          sparse_x_dense(j1,:) = sparse_x_dense(j1,:) + sparse(j1,j2)*dense(j2,:)
     
    376376      mblock = m/3
    377377      m2block = 2*mblock
    378       do j4 = 1,nrepeat
     378      DO j4 = 1,nrepeat
    379379        repeated_square = 0.0_jprb
    380380        ! Do the top-left (C, D, F & G)
    381         do j2 = 1,m2block
    382           do j1 = 1,m2block
    383             do j3 = 1,m2block
     381        DO j2 = 1,m2block
     382          DO j1 = 1,m2block
     383            DO j3 = 1,m2block
    384384              repeated_square(j1,j2) = repeated_square(j1,j2) &
    385385                   &                 + A(j1,j3)*A(j3,j2)
     
    387387          end do
    388388        end do
    389         do j2 = m2block+1, m
     389        DO j2 = m2block+1, m
    390390          ! Do the top-right (E & H)
    391           do j1 = 1,m2block
    392             do j3 = 1,m
     391          DO j1 = 1,m2block
     392            DO j3 = 1,m
    393393              repeated_square(j1,j2) = repeated_square(j1,j2) &
    394394                   &                 + A(j1,j3)*A(j3,j2)
     
    396396          end do
    397397          ! Do the bottom-right (I)
    398           do j1 = m2block+1, m
    399             do j3 = m2block+1,m
     398          DO j1 = m2block+1, m
     399            DO j3 = m2block+1,m
    400400              repeated_square(j1,j2) = repeated_square(j1,j2) &
    401401                   &                 + A(j1,j3)*A(j3,j2)
     
    409409    else
    410410      ! Ordinary dense matrix
    411       do j4 = 1,nrepeat
     411      DO j4 = 1,nrepeat
    412412        repeated_square = 0.0_jprb
    413         do j2 = 1,m
    414           do j1 = 1,m
    415             do j3 = 1,m
     413        DO j2 = 1,m
     414          DO j1 = 1,m
     415            DO j3 = 1,m
    416416              repeated_square(j1,j2) = repeated_square(j1,j2) &
    417417                   &                 + A(j1,j3)*A(j3,j2)
     
    549549    U33 = A(1:iend,3,3) - L31*A(1:iend,1,3) - L32*U23
    550550
    551     do j = 1,3
     551    DO j = 1,3
    552552      ! Solve Ly = B(:,:,j) by forward substitution
    553553      ! y1 = B(:,1,j)
     
    649649    LU(1:iend,1:m,1:m) = A(1:iend,1:m,1:m)
    650650
    651     do j2 = 1, m
    652       do j1 = 1, j2-1
     651    DO j2 = 1, m
     652      DO j1 = 1, j2-1
    653653        s = LU(1:iend,j1,j2)
    654         do j3 = 1, j1-1
     654        DO j3 = 1, j1-1
    655655          s = s - LU(1:iend,j1,j3) * LU(1:iend,j3,j2)
    656656        end do
    657657        LU(1:iend,j1,j2) = s
    658658      end do
    659       do j1 = j2, m
     659      DO j1 = j2, m
    660660        s = LU(1:iend,j1,j2)
    661         do j3 = 1, j2-1
     661        DO j3 = 1, j2-1
    662662          s = s - LU(1:iend,j1,j3) * LU(1:iend,j3,j2)
    663663        end do
     
    666666      if (j2 /= m) then
    667667        s = 1.0_jprb / LU(1:iend,j2,j2)
    668         do j1 = j2+1, m
     668        DO j1 = j2+1, m
    669669          LU(1:iend,j1,j2) = LU(1:iend,j1,j2) * s
    670670        end do
     
    691691
    692692    ! First solve Ly=b
    693     do j2 = 2, m
    694       do j1 = 1, j2-1
     693    DO j2 = 2, m
     694      DO j1 = 1, j2-1
    695695        x(1:iend,j2) = x(1:iend,j2) - x(1:iend,j1)*LU(1:iend,j2,j1)
    696696      end do
    697697    end do
    698698    ! Now solve Ux=y
    699     do j2 = m, 1, -1
    700       do j1 = j2+1, m
     699    DO j2 = m, 1, -1
     700      DO j1 = j2+1, m
    701701        x(1:iend,j2) = x(1:iend,j2) - x(1:iend,j1)*LU(1:iend,j2,j1)
    702702      end do
     
    722722    call lu_factorization(n,iend,m,A,LU)
    723723
    724     do j = 1, m
     724    DO j = 1, m
    725725      call lu_substitution(n,iend,m,LU,B(1:,1:m,j),X(1:iend,1:m,j))
    726726!      call lu_substitution(n,iend,m,LU,B(1:n,1:m,j),X(1:iend,1:m,j))
     
    835835
    836836    ! Compute the 1-norms of A
    837     do j3 = 1,m
     837    DO j3 = 1,m
    838838      sum_column(:) = 0.0_jprb
    839       do j2 = 1,m
    840         do j1 = 1,iend
     839      DO j2 = 1,m
     840        DO j1 = 1,iend
    841841          sum_column(j1) = sum_column(j1) + abs(A(j1,j2,j3))
    842842        end do
    843843      end do
    844       do j1 = 1,iend
     844      DO j1 = 1,iend
    845845        if (sum_column(j1) > normA(j1)) then
    846846          normA(j1) = sum_column(j1)
     
    861861    ! Scale the input matrices by a power of 2
    862862    scaling = 2.0_jprb**(-expo)
    863     do j3 = 1,m
    864       do j2 = 1,m
     863    DO j3 = 1,m
     864      DO j2 = 1,m
    865865        A(1:iend,j2,j3) = A(1:iend,j2,j3) * scaling
    866866      end do
     
    872872
    873873    V = c(8)*A6 + c(6)*A4 + c(4)*A2
    874     do j3 = 1,m
     874    DO j3 = 1,m
    875875      V(:,j3,j3) = V(:,j3,j3) + c(2)
    876876    end do
     
    878878    V = c(7)*A6 + c(5)*A4 + c(3)*A2
    879879    ! Add a multiple of the identity matrix
    880     do j3 = 1,m
     880    DO j3 = 1,m
    881881      V(:,j3,j3) = V(:,j3,j3) + c(1)
    882882    end do
     
    887887
    888888    ! Add the identity matrix
    889     do j3 = 1,m
     889    DO j3 = 1,m
    890890      A(1:iend,j3,j3) = A(1:iend,j3,j3) + 1.0_jprb
    891891    end do
    892892
    893893    ! Loop through the matrices
    894     do j1 = 1,iend
     894    DO j1 = 1,iend
    895895      if (expo(j1) > 0) then
    896896        ! Square matrix j1 expo(j1) times         
     
    10161016
    10171017    ! Compute V * diag_rdivide_V
    1018     do j1 = 1,3
    1019       do j2 = 1,3
     1018    DO j1 = 1,3
     1019      DO j2 = 1,3
    10201020        R(1:iend,j2,j1) = V(1:iend,j2,1)*diag_rdivide_V(1:iend,1,j1) &
    10211021             &          + V(1:iend,j2,2)*diag_rdivide_V(1:iend,2,j1) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_mcica_lw.F90

    r4946 r5158  
    146146
    147147    ! Loop through columns
    148     do jcol = istartcol,iendcol
     148    DO jcol = istartcol,iendcol
    149149
    150150      ! Clear-sky calculation
     
    203203        is_clear_sky_layer = .true.
    204204        i_cloud_top = nlev+1
    205         do jlev = 1,nlev
     205        DO jlev = 1,nlev
    206206          ! Compute combined gas+aerosol+cloud optical properties
    207207          if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
     
    212212            end if
    213213
    214             do jg = 1,ng
     214            DO jg = 1,ng
    215215              od_cloud_new(jg) = od_scaling(jg,jlev) &
    216216                 &  * od_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol)
     
    228228                ! case that od_total > 0.0 and ssa_total > 0.0 but
    229229                ! od_total*ssa_total == 0 due to underflow
    230                 do jg = 1,ng
     230                DO jg = 1,ng
    231231                  if (od_total(jg) > 0.0_jprb) then
    232232                    scat_od_total(jg) = ssa(jg,jlev,jcol)*od(jg,jlev,jcol) &
     
    247247              else
    248248
    249                 do jg = 1,ng
     249                DO jg = 1,ng
    250250                  if (od_total(jg) > 0.0_jprb) then
    251251                    scat_od = ssa_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90

    r4946 r5158  
    138138
    139139    ! Loop through columns
    140     do jcol = istartcol,iendcol
     140    DO jcol = istartcol,iendcol
    141141      ! Only perform calculation if sun above the horizon
    142142      if (single_level%cos_sza(jcol) > 0.0_jprb) then
     
    155155        else
    156156          ! Apply delta-Eddington scaling to the aerosol-gas mixture
    157           do jlev = 1,nlev
     157          DO jlev = 1,nlev
    158158            od_total  =  od(:,jlev,jcol)
    159159            ssa_total = ssa(:,jlev,jcol)
     
    205205        if (total_cloud_cover >= config%cloud_fraction_threshold) then
    206206          ! Total-sky calculation
    207           do jlev = 1,nlev
     207          DO jlev = 1,nlev
    208208            ! Compute combined gas+aerosol+cloud optical properties
    209209            if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
    210               do jg = 1,ng
     210              DO jg = 1,ng
    211211                od_cloud_new(jg) = od_scaling(jg,jlev) &
    212212                   &  * od_cloud(config%i_band_from_reordered_g_sw(jg),jlev,jcol)
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_monochromatic.F90

    r4773 r5158  
    164164    integer :: jlev
    165165
    166     do jlev = 1,nlev
     166    DO jlev = 1,nlev
    167167      ! The fraction of the total optical depth in the current layer
    168168      ! is proportional to the fraction of the mass of the atmosphere
     
    192192             &  = StefanBoltzmann * single_level%skin_temperature(istartcol:iendcol)**4 &
    193193             &  * single_level%lw_emissivity(istartcol:iendcol,1)
    194         do jlev = 1,nlev+1
     194        DO jlev = 1,nlev+1
    195195          planck_hl(1,jlev,istartcol:iendcol) = StefanBoltzmann * thermodynamics%temperature_hl(istartcol:iendcol,jlev)**4
    196196        end do
     
    200200             &             single_level%skin_temperature(istartcol:iendcol)) &
    201201             &  * single_level%lw_emissivity(istartcol:iendcol,1)
    202         do jlev = 1,nlev+1
     202        DO jlev = 1,nlev+1
    203203          planck_hl(1,jlev,istartcol:iendcol) = Pi*planck_function(config%mono_lw_wavelength, &
    204204               &             thermodynamics%temperature_hl(istartcol:iendcol,jlev))
     
    261261    ! Convert cloud mixing ratio into liquid and ice water path
    262262    ! in each layer
    263     do jlev = 1, nlev
    264       do jcol = istartcol, iendcol
     263    DO jlev = 1, nlev
     264      DO jcol = istartcol, iendcol
    265265        ! Factor to convert from gridbox-mean mass mixing ratio to
    266266        ! in-cloud water path involves the pressure difference in
     
    288288
    289289    if (config%iverbose >= 4) then
    290       do jcol = istartcol,iendcol
     290      DO jcol = istartcol,iendcol
    291291        write(*,'(a,i0,a,f7.3,a,f7.3)') 'Profile ', jcol, ': shortwave optical depth = ', &
    292292             &  sum(od_sw_cloud(1,:,jcol)*cloud%fraction(jcol,:)), &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_overlap.F90

    r4773 r5158  
    9494    ! multiply by "factor"
    9595    denominator = 1.0_jprb
    96     do jreg = 1,nreg
     96    DO jreg = 1,nreg
    9797      op_x_frac_min(jreg) = op(jreg) &
    9898           &  * min(frac_upper(jreg), frac_lower(jreg))
     
    103103      factor = 1.0_jprb / denominator
    104104      ! Create the random part of the overlap matrix
    105       do jupper = 1,nreg
    106         do jlower = 1,nreg
     105      DO jupper = 1,nreg
     106        DO jlower = 1,nreg
    107107          overlap_matrix(jupper,jlower) = factor &
    108108               &  * (frac_lower(jlower)-op_x_frac_min(jlower)) &
     
    115115   
    116116    ! Add on the maximum part of the overlap matrix
    117     do jreg = 1,nreg
     117    DO jreg = 1,nreg
    118118      overlap_matrix(jreg,jreg) = overlap_matrix(jreg,jreg) &
    119119           &  + op_x_frac_min(jreg)
     
    368368
    369369    ! Loop through each atmospheric column
    370     do jcol = istartcol, iendcol
     370    DO jcol = istartcol, iendcol
    371371      ! For this column, outer space is treated as one clear-sky
    372372      ! region, so the fractions are assigned as such
     
    381381      ! half-level starting at 1 for the top-of-atmosphere, as well
    382382      ! as indexing each level starting at 1 for the top-most level.
    383       do jlev = 1,nlev+1
     383      DO jlev = 1,nlev+1
    384384        ! Fraction of each region just below the interface
    385385        if (jlev > nlev) then
     
    426426
    427427        ! Convert to directional overlap matrices
    428         do jupper = 1,nreg
    429           do jlower = 1,nreg
     428        DO jupper = 1,nreg
     429          DO jlower = 1,nreg
    430430            if (frac_lower(jlower) >= frac_threshold) then
    431431              u_matrix(jupper,jlower,jlev,jcol) = overlap_matrix(jupper,jlower) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_pdf_sampler.F90

    r4773 r5158  
    189189    real(jprb) :: wfsd, wcdf
    190190
    191     do jsamp = 1,nsamp
     191    DO jsamp = 1,nsamp
    192192      if (mask(jsamp)) then
    193193        ! Bilinear interpolation with bounds
     
    237237    real(jprb) :: wfsd, wcdf
    238238
    239     do jz = 1,nz
    240       do jg = 1,ng
     239    DO jz = 1,nz
     240      DO jg = 1,ng
    241241        if (cdf(jg, jz) > 0.0_jprb) then
    242242          ! Bilinear interpolation with bounds
     
    291291    real(jprb) :: wfsd, wcdf
    292292
    293     do jz = 1,nz
     293    DO jz = 1,nz
    294294
    295295      if (mask(jz)) then
    296296       
    297         do jg = 1,ng
     297        DO jg = 1,ng
    298298          if (cdf(jg, jz) > 0.0_jprb) then
    299299            ! Bilinear interpolation with bounds
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_random_numbers.F90

    r4773 r5158  
    165165      ! populate the state
    166166      rseed = real(abs(this%iseed),jprd)
    167       do jstr = 1,this%nmaxstreams
     167      DO jstr = 1,this%nmaxstreams
    168168        ! Note that nint returns an integer of type jpib (8-byte)
    169169        ! which may be converted to double if that is the type of
     
    174174
    175175      ! One warmup of the C++ minstd_rand algorithm
    176       do jstr = 1,this%nmaxstreams
     176      DO jstr = 1,this%nmaxstreams
    177177        this%istate(jstr) = mod(IMinstdA * this%istate(jstr), IMinstdM)
    178178      end do
     
    182182      call random_seed(size=nseed)
    183183      allocate(iseednative(nseed))
    184       do jseed = 1,nseed
     184      DO jseed = 1,nseed
    185185        iseednative(jseed) = this%iseed + jseed - 1
    186186      end do
     
    208208
    209209      ! C++ minstd_rand algorithm
    210       do i = 1, imax
     210      DO i = 1, imax
    211211        ! The following calculation is computed entirely with 8-byte
    212212        ! numbers (whether real or integer)
     
    243243
    244244      ! C++ minstd_ran algorithm
    245       do jblock = 1,size(randnum,2)
     245      DO jblock = 1,size(randnum,2)
    246246        ! These lines should be vectorizable
    247         do i = 1, imax
     247        DO i = 1, imax
    248248          this%istate(i) = mod(IMinstdA * this%istate(i), IMinstdM)
    249249          randnum(i,jblock) = IMinstdScale * this%istate(i)
     
    278278
    279279      ! C++ minstd_ran algorithm
    280       do jblock = 1,size(randnum,2)
     280      DO jblock = 1,size(randnum,2)
    281281        if (mask(jblock)) then
    282282          ! These lines should be vectorizable
    283           do i = 1, imax
     283          DO i = 1, imax
    284284            this%istate(i) = mod(IMinstdA * this%istate(i), IMinstdM)
    285285            randnum(i,jblock) = IMinstdScale * this%istate(i)
     
    290290    else
    291291
    292       do jblock = 1,size(randnum,2)
     292      DO jblock = 1,size(randnum,2)
    293293        if (mask(jblock)) then
    294294          call random_number(randnum(:,jblock))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_regions.F90

    r4773 r5158  
    122122        ! sigma, then the 16th percentile of the lognormal is very
    123123        ! close to exp(mu-sigma).
    124         do jcol = istartcol,iendcol
    125           do jlev = 1,nlev
     124        DO jcol = istartcol,iendcol
     125          DO jlev = 1,nlev
    126126            if (cloud_fraction(jcol,jlev) < frac_threshold) then
    127127              reg_fracs(1,jlev,jcol)       = 1.0_jprb
     
    145145        ! cover of the denser region - see region_fractions routine
    146146        ! below
    147         do jcol = istartcol,iendcol
    148           do jlev = 1,nlev
     147        DO jcol = istartcol,iendcol
     148          DO jlev = 1,nlev
    149149            if (cloud_fraction(jcol,jlev) < frac_threshold) then
    150150              reg_fracs(1,jlev,jcol)    = 1.0_jprb
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_save.F90

    r4853 r5158  
    11681168         &   dim2_name="column", dim1_name="level", &
    11691169         &   units_str="1", long_name="Ozone mass mixing ratio")
    1170     do jgas = 1,NMaxGases
     1170    DO jgas = 1,NMaxGases
    11711171      if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
    11721172        write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
     
    12591259    call gas%get(IO3, IMassMixingRatio, mixing_ratio)
    12601260    call out_file%put("o3_mmr", mixing_ratio)
    1261     do jgas = 1,NMaxGases
     1261    DO jgas = 1,NMaxGases
    12621262      if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
    12631263        write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_single_level.F90

    r4773 r5158  
    206206    end if
    207207
    208     do jcol = istartcol,iendcol
     208    DO jcol = istartcol,iendcol
    209209      this%iseed(jcol) = jcol
    210210    end do
     
    275275
    276276        sw_albedo_band = 0.0_jprb
    277         do jband = 1,config%n_bands_sw
    278           do jalbedoband = 1,nalbedoband
     277        DO jband = 1,config%n_bands_sw
     278          DO jalbedoband = 1,nalbedoband
    279279            if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
    280               do jcol = istartcol,iendcol
     280              DO jcol = istartcol,iendcol
    281281                sw_albedo_band(jcol,jband) &
    282282                    &  = sw_albedo_band(jcol,jband) &
     
    292292        if (allocated(this%sw_albedo_direct)) then
    293293          sw_albedo_band = 0.0_jprb
    294           do jband = 1,config%n_bands_sw
    295             do jalbedoband = 1,nalbedoband
     294          DO jband = 1,config%n_bands_sw
     295            DO jalbedoband = 1,nalbedoband
    296296              if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
    297297                sw_albedo_band(istartcol:iendcol,jband) &
     
    342342        end if
    343343        lw_albedo_band = 0.0_jprb
    344         do jband = 1,config%n_bands_lw
    345           do jalbedoband = 1,nalbedoband
     344        DO jband = 1,config%n_bands_lw
     345          DO jalbedoband = 1,nalbedoband
    346346            if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then
    347               do jcol = istartcol,iendcol
     347              DO jcol = istartcol,iendcol
    348348                lw_albedo_band(jcol,jband) &
    349349                    &  = lw_albedo_band(jcol,jband) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_lw.F90

    r4773 r5158  
    313313
    314314    ! Main loop over columns
    315     do jcol = istartcol, iendcol
     315    DO jcol = istartcol, iendcol
    316316      ! --------------------------------------------------------
    317317      ! Section 2: Prepare column-specific variables and arrays
     
    325325      ! cloud%crop_cloud_fraction has already been called
    326326      is_clear_sky_layer = .true.
    327       do jlev = 1,nlev
     327      DO jlev = 1,nlev
    328328        if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
    329329          is_clear_sky_layer(jlev) = .false.
     
    336336      ! In this section the reflectance, transmittance and sources
    337337      ! are computed for each layer
    338       do jlev = 1,nlev ! Start at top-of-atmosphere
     338      DO jlev = 1,nlev ! Start at top-of-atmosphere
    339339        ! --------------------------------------------------------
    340340        ! Section 3.1: Layer-specific initialization
     
    395395            ! approximate order of gas optical depth.
    396396            ng3D = ng
    397             do jg = 1, ng
     397            DO jg = 1, ng
    398398              if (od_region(jg,1) > config%max_gas_od_3D) then
    399399                ng3D = jg-1
     
    488488              end if
    489489
    490               do jreg = 1, nreg-1
     490              DO jreg = 1, nreg-1
    491491                ! Compute lateral transfer rate from region jreg to
    492492                ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but
     
    536536          ! g-point, mapping from the cloud properties
    537537          ! defined in each band.
    538           do jg = 1,ng
     538          DO jg = 1,ng
    539539            ! Mapping from g-point to band
    540540            iband = config%i_band_from_reordered_g_lw(jg)
     
    544544
    545545            ! Loop over cloudy regions
    546             do jreg = 2,nreg
     546            DO jreg = 2,nreg
    547547              ! Add scaled cloud optical depth to clear-sky value
    548548              od_region(jg,jreg) = od_region(jg,1) &
     
    603603          ! exponential, then computing the various
    604604          ! transmission/reflectance matrices from that.
    605           do jreg = 1,nregActive
     605          DO jreg = 1,nregActive
    606606            ! Write the diagonal elements of -Gamma1*z1
    607607            Gamma_z1(1:ng3D,jreg,jreg) &
     
    628628            ! non-zeros in the cloudy-sky regions even if the cloud
    629629            ! fraction is zero
    630             do jreg = nregActive+1,nreg
     630            DO jreg = nregActive+1,nreg
    631631              Gamma_z1(1:ng3D,jreg,jreg) = Gamma_z1(1:ng3D,1,1)
    632632              Gamma_z1(1:ng3D,nreg+jreg,jreg) = Gamma_z1(1:ng3D,nreg+1,1)
     
    656656          end if
    657657
    658           do jreg = 1,nregActive-1
     658          DO jreg = 1,nregActive-1
    659659            ! Add the terms assocated with 3D transport to Gamma1*z1.
    660660            ! First the rate from region jreg to jreg+1
     
    775775          ! formulas for reflectance and transmittance, and equivalent
    776776          ! solutions to the coupled ODEs for the sources.
    777           do jreg = 2, nregActive
     777          DO jreg = 2, nregActive
    778778            call calc_reflectance_transmittance_lw(ng-ng3D, &
    779779                 &  od_region(ng3D+1:ng,jreg), &
     
    806806      ! Calculate the upwelling radiation emitted by the surface, and
    807807      ! copy the surface albedo into total_albedo
    808       do jreg = 1,nreg
    809         do jg = 1,ng
     808      DO jreg = 1,nreg
     809        DO jg = 1,ng
    810810          ! region_fracs(jreg,nlev,jcol) is the fraction of each
    811811          ! region in the lowest model level
     
    816816      ! Equivalent surface values for computing clear-sky fluxes
    817817      if (config%do_clear) then
    818         do jg = 1,ng
     818        DO jg = 1,ng
    819819          total_source_clear(jg,nlev+1) = emission(jg,jcol)
    820820        end do
     
    826826      ! Loop back up through the atmosphere computing the total albedo
    827827      ! and the total upwelling due to emission below each level
    828       do jlev = nlev,1,-1
     828      DO jlev = nlev,1,-1
    829829        if (config%do_clear) then
    830830          ! Use adding method for clear-sky arrays; note that there
     
    876876          total_albedo_below = 0.0_jprb
    877877          total_source_below = 0.0_jprb
    878           do jreg = 1,nreg
     878          DO jreg = 1,nreg
    879879            inv_denom_scalar(:) = 1.0_jprb / (1.0_jprb &
    880880                 &  - total_albedo(:,jreg,jreg,jlev+1)*reflectance(:,jreg,jreg,jlev))
     
    918918            ! essentially diag(total_albedo) = matmul(transpose(v_matrix),
    919919            ! diag(total_albedo_below)).
    920             do jreg = 1,nreg
    921               do jreg2 = 1,nreg
     920            DO jreg = 1,nreg
     921              DO jreg2 = 1,nreg
    922922                total_albedo(:,jreg,jreg,jlev) &
    923923                     &  = total_albedo(:,jreg,jreg,jlev) &
     
    966966
    967967      ! Final loop back down through the atmosphere to compute fluxes
    968       do jlev = 1,nlev
     968      DO jlev = 1,nlev
    969969        if (config%do_clear) then
    970970          ! Scalar operations for clear-sky fluxes
     
    999999               &  flux_dn_above) + total_source(:,:,jlev+1)
    10001000        else
    1001           do jreg = 1,nreg
     1001          DO jreg = 1,nreg
    10021002            ! Scalar operations for all regions, requiring that
    10031003            ! reflectance, transmittance and total_albedo are diagonal
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_sw.F90

    r4853 r5158  
    328328
    329329    ! Main loop over columns
    330     do jcol = istartcol, iendcol
     330    DO jcol = istartcol, iendcol
    331331      ! --------------------------------------------------------
    332332      ! Section 2: Prepare column-specific variables and arrays
     
    404404      is_clear_sky_layer = .true.
    405405      i_cloud_top = nlev+1
    406       do jlev = nlev,1,-1
     406      DO jlev = nlev,1,-1
    407407        if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
    408408          is_clear_sky_layer(jlev) = .false.
     
    416416      ! In this section the reflectance, transmittance and sources
    417417      ! are computed for each layer
    418       do jlev = 1,nlev ! Start at top-of-atmosphere
     418      DO jlev = 1,nlev ! Start at top-of-atmosphere
    419419        ! --------------------------------------------------------
    420420        ! Section 3.1: Layer-specific initialization
     
    469469            ! approximate order of gas optical depth.
    470470            ng3D = ng
    471             do jg = 1, ng
     471            DO jg = 1, ng
    472472              if (od_region(jg,1) > config%max_gas_od_3D) then
    473473                ng3D = jg-1
     
    553553              end if
    554554
    555               do jreg = 1, nreg-1
     555              DO jreg = 1, nreg-1
    556556                ! Compute lateral transfer rates from region jreg to
    557557                ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but
     
    615615          ! g-point, mapping from the cloud properties
    616616          ! defined in each band.
    617           do jg = 1,ng
     617          DO jg = 1,ng
    618618            ! Mapping from g-point to band
    619619            iband = config%i_band_from_reordered_g_sw(jg)
     
    629629
    630630            ! Loop over cloudy regions
    631             do jreg = 2,nreg
     631            DO jreg = 2,nreg
    632632              scat_od_cloud = od_cloud(iband,jlev,jcol) &
    633633                   &  * ssa_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
     
    679679          ! exponential, then computing the various
    680680          ! transmission/reflectance matrices from that.
    681           do jreg = 1,nregactive
     681          DO jreg = 1,nregactive
    682682            ! Write the diagonal elements of -Gamma1*z1
    683683            Gamma_z1(1:ng3D,jreg,jreg) &
     
    701701          end do
    702702
    703           do jreg = 1,nregactive-1
     703          DO jreg = 1,nregactive-1
    704704            ! Write the elements of -Gamma1*z1 concerned with 3D
    705705            ! transport
     
    820820          ! Compute reflectance, transmittance and associated terms
    821821          ! for each cloudy region, using the Meador-Weaver formulas
    822           do jreg = 2, nregactive
     822          DO jreg = 2, nregactive
    823823            call calc_reflectance_transmittance_sw(ng-ng3D, &
    824824                 &  mu0, &
     
    853853      ! beam incident on the surface, and copy the surface albedo
    854854      ! into total_albedo
    855       do jreg = 1,nreg
    856         do jg = 1,ng
     855      DO jreg = 1,nreg
     856        DO jg = 1,ng
    857857          total_albedo(jg,jreg,jreg,nlev+1) = albedo_diffuse(jg,jcol)
    858858          total_albedo_direct(jg,jreg,jreg,nlev+1) &
     
    875875      ! Work back up through the atmosphere computing the total albedo
    876876      ! of the atmosphere below that point using the adding method
    877       do jlev = nlev,1,-1
     877      DO jlev = nlev,1,-1
    878878
    879879        ! --------------------------------------------------------
     
    995995          ! radiation:
    996996          total_albedo(:,:,:,jlev) = 0.0_jprb
    997           do jreg = 1,nreg    ! Target layer (jlev-1)
    998             do jreg2 = 1,nreg ! Current layer (jlev)
     997          DO jreg = 1,nreg    ! Target layer (jlev-1)
     998            DO jreg2 = 1,nreg ! Current layer (jlev)
    999999              total_albedo(:,jreg,jreg,jlev) = total_albedo(:,jreg,jreg,jlev) &
    10001000                   &  + sum(total_albedo_below(:,:,jreg2),2) &
     
    10041004          ! ...then direct radiation:
    10051005          total_albedo_direct(:,:,:,jlev) = 0.0_jprb
    1006           do jreg = 1,nreg    ! Target layer (jlev-1)
    1007             do jreg2 = 1,nreg ! Current layer (jlev)
     1006          DO jreg = 1,nreg    ! Target layer (jlev-1)
     1007            DO jreg2 = 1,nreg ! Current layer (jlev)
    10081008              total_albedo_direct(:,jreg,jreg,jlev) = total_albedo_direct(:,jreg,jreg,jlev) &
    10091009                   &  + sum(total_albedo_below_direct(:,:,jreg2),2) &
     
    10311031          ! First diffuse radiation:
    10321032          albedo_part = total_albedo_below
    1033           do jreg = 1,nreg
     1033          DO jreg = 1,nreg
    10341034            albedo_part(:,jreg,jreg) = 0.0_jprb
    10351035          end do
     
    10401040          ! ...then direct radiation:
    10411041          albedo_part = total_albedo_below_direct
    1042           do jreg = 1,nreg
     1042          DO jreg = 1,nreg
    10431043            albedo_part(:,jreg,jreg) = 0.0_jprb
    10441044          end do
     
    10591059            ! essentially diag(total_albedo) += matmul(transpose(v_matrix),
    10601060            ! diag(total_albedo_below)).
    1061             do jreg = 1,nreg
    1062               do jreg2 = 1,nreg
     1061            DO jreg = 1,nreg
     1062              DO jreg2 = 1,nreg
    10631063                total_albedo(:,jreg,jreg,jlev) &
    10641064                     &  = total_albedo(:,jreg,jreg,jlev) &
     
    10751075            ! "Explicit entrapment"
    10761076
    1077             do jreg2 = 1,nreg
     1077            DO jreg2 = 1,nreg
    10781078              ! Loop through each region in the lower layer. For one
    10791079              ! of the regions in the lower layer, we are imagining it
     
    11101110                     &  / max(config%cloud_fraction_threshold, region_fracs(jreg2,jlev,jcol))
    11111111
    1112                 do jreg = 1, nreg-1
     1112                DO jreg = 1, nreg-1
    11131113                  ! Compute lateral transfer rates from region jreg to
    11141114                  ! jreg+1 as before, but without length scale which
     
    11391139              inv_effective_size = min(cloud%inv_cloud_effective_size(jcol,jlev-1), &
    11401140                   &                   1.0_jprb/config%min_cloud_effective_size)
    1141               do jreg = 1,nreg-1
     1141              DO jreg = 1,nreg-1
    11421142                ! Diffuse transport down and up with one random
    11431143                ! scattering event
     
    11671167              ! trival limit, so we cap the maximum input to expm by
    11681168              ! scaling down if necessary
    1169               do jg = 1,ng
     1169              DO jg = 1,ng
    11701170                max_entr = -min(entrapment(jg,1,1),entrapment(jg,2,2))
    11711171                if (max_entr > config%max_cloud_od) then
     
    11971197#ifndef EXPLICIT_EDGE_ENTRAPMENT
    11981198              ! Scale to get the contribution to the diffuse albedo
    1199               do jreg3 = 1,nreg
    1200                 do jreg = 1,nreg
     1199              DO jreg3 = 1,nreg
     1200                DO jreg = 1,nreg
    12011201                  albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) &
    12021202                       &  * v_matrix(jreg2,jreg,jlev,jcol) * total_albedo_below(:,jreg2,jreg2)
     
    12101210              albedo_part = 0.0_jprb
    12111211              ! Scale to get the contribution to the diffuse albedo
    1212               do jreg3 = 1,nreg     ! TO upper region
    1213                 do jreg = 1,nreg    ! FROM upper region
     1212              DO jreg3 = 1,nreg     ! TO upper region
     1213                DO jreg = 1,nreg    ! FROM upper region
    12141214                  transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) &
    12151215                       &  * cloud%overlap_param(jcol,jlev-1) &
    12161216                       &  * min(region_fracs(jreg,jlev,jcol), region_fracs(jreg,jlev,jcol)) &
    12171217                       &  / max(config%cloud_fraction_threshold, region_fracs(jreg,jlev,jcol))
    1218                   do jreg4 = 1,nreg ! VIA first lower region (jreg2 is second lower region)
     1218                  DO jreg4 = 1,nreg ! VIA first lower region (jreg2 is second lower region)
    12191219                    if (.not. (jreg4 == jreg .and. jreg4 /= jreg2)) then
    12201220                      albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) + entrapment(:,jreg3,jreg) &
     
    12361236              ! Now do the same for the direct albedo
    12371237              entrapment = 0.0_jprb
    1238               do jreg = 1,nreg-1
     1238              DO jreg = 1,nreg-1
    12391239                ! Direct transport down and diffuse up with one random
    12401240                ! scattering event
     
    12641264              ! trival limit, so we cap the maximum input to expm by
    12651265              ! scaling down if necessary
    1266               do jg = 1,ng
     1266              DO jg = 1,ng
    12671267                max_entr = -min(entrapment(jg,1,1),entrapment(jg,2,2))
    12681268                if (max_entr > config%max_cloud_od) then
     
    12891289
    12901290#ifndef EXPLICIT_EDGE_ENTRAPMENT
    1291               do jreg3 = 1,nreg
    1292                 do jreg = 1,nreg
     1291              DO jreg3 = 1,nreg
     1292                DO jreg = 1,nreg
    12931293                  albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) &
    12941294                       &  * v_matrix(jreg2,jreg,jlev,jcol) * total_albedo_below_direct(:,jreg2,jreg2)
     
    12981298              entrapment = albedo_part
    12991299              albedo_part = 0.0_jprb
    1300               do jreg3 = 1,nreg
    1301                 do jreg = 1,nreg
     1300              DO jreg3 = 1,nreg
     1301                DO jreg = 1,nreg
    13021302                  transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) &
    13031303                       &  * cloud%overlap_param(jcol,jlev-1) &
    13041304                       &  * min(region_fracs(jreg,jlev,jcol), region_fracs(jreg,jlev-1,jcol)) &
    13051305                       &  / max(config%cloud_fraction_threshold, region_fracs(jreg,jlev,jcol))
    1306                   do jreg4 = 1,nreg
     1306                  DO jreg4 = 1,nreg
    13071307                    if (.not. (jreg4 == jreg .and. jreg4 /= jreg2)) then
    13081308                     albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) + entrapment(:,jreg3,jreg) &
     
    13431343          end if
    13441344
    1345           do jreg = 1,nreg          ! Target layer (jlev-1)
    1346             do jreg2 = 1,nregactive ! Current layer (jlev)
     1345          DO jreg = 1,nreg          ! Target layer (jlev-1)
     1346            DO jreg2 = 1,nregactive ! Current layer (jlev)
    13471347              x_direct_above(:,jreg) = x_direct_above(:,jreg) &
    13481348                   &  + x_direct(:,jreg2) * v_matrix(jreg2,jreg,jlev,jcol)
     
    13681368      ! Direct downwelling flux (into a plane perpendicular to the
    13691369      ! sun) entering the top of each region in the top-most layer
    1370       do jreg = 1,nreg
     1370      DO jreg = 1,nreg
    13711371        direct_dn_below(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol)
    13721372      end do
     
    14211421
    14221422      ! Final loop back down through the atmosphere to compute fluxes
    1423       do jlev = 1,nlev
     1423      DO jlev = 1,nlev
    14241424
    14251425#ifdef PRINT_ENTRAPMENT_DATA
     
    16751675         &                             + tan_diffuse_angle_3d*tan_diffuse_angle_3d) * 0.5_jprb
    16761676
    1677     do jreg = istartreg,iendreg
     1677    DO jreg = istartreg,iendreg
    16781678      ! Geometric series enhancement due to multiple scattering: the
    16791679      ! amplitude enhancement is equal to the limit of
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spectral_definition.F90

    r4853 r5158  
    201201    else
    202202      find_wavenumber = 1
    203       do while (wavenumber > this%wavenumber2(find_wavenumber) &
     203      DO while (wavenumber > this%wavenumber2(find_wavenumber) &
    204204           &    .and. find_wavenumber < this%nwav)
    205205        find_wavenumber = find_wavenumber + 1
     
    284284      end if
    285285
    286       do jband = 1,this%nband
     286      DO jband = 1,this%nband
    287287        weight = 0.0_jprb
    288         do jwav = 1,nwav
     288        DO jwav = 1,nwav
    289289          ! Work out wavenumber range for which this cloud wavenumber
    290290          ! will be applicable
     
    321321            ! Find interpolating points
    322322            iwav = 2
    323             do while (wavenumber(iwav) < this%wavenumber2_band(jband))
     323            DO while (wavenumber(iwav) < this%wavenumber2_band(jband))
    324324              iwav = iwav+1
    325325            end do
     
    358358      mapping = 0.0_jprb
    359359      ! Loop over wavenumbers representing cloud
    360       do jwav = 1,nwav
     360      DO jwav = 1,nwav
    361361        ! Clear the weights. The weight says for one wavenumber in the
    362362        ! cloud file, what is its fractional contribution to each of
     
    400400                 &  / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
    401401            if (isd1-isd0 > 1) then
    402               do isd = isd0+1,isd1-1
     402              DO isd = isd0+1,isd1-1
    403403                ! Intermediate trapezia
    404404                weight(isd) = 0.5_jprb * (this%wavenumber1(isd)+this%wavenumber2(isd) &
     
    447447                 &  / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
    448448            if (isd2-isd1 > 1) then
    449               do isd = isd1+1,isd2-1
     449              DO isd = isd1+1,isd2-1
    450450                ! Intermediate trapezia
    451451                weight(isd) = weight(isd) + 0.5_jprb * (2.0_jprb*wavenum2 &
     
    468468        weight = weight * planck_weight
    469469
    470         do jg = 1,this%ng
     470        DO jg = 1,this%ng
    471471          mapping(jg, jwav) = sum(weight * this%gpoint_fraction(:,jg))
    472472        end do
     
    478478
    479479      ! Normalize mapping matrix
    480       do jg = 1,this%ng
     480      DO jg = 1,this%ng
    481481        mapping(jg,:) = mapping(jg,:) * (1.0_jprb/sum(mapping(jg,:)))
    482482      end do
     
    589589    ! Check wavelength is monotonically increasing
    590590    if (ninterval > 2) then
    591       do jint = 2,ninterval-1
     591      DO jint = 2,ninterval-1
    592592        if (wavelength_bound(jint) <= wavelength_bound(jint-1)) then
    593593          write(nulerr, '(a)') '*** Error: wavelength bounds must be monotonically increasing'
     
    609609      end if
    610610
    611       do jband = 1,this%nband
    612         do jint = 1,ninterval
     611      DO jband = 1,this%nband
     612        DO jint = 1,ninterval
    613613          if (jint == 1) then
    614614            ! First input interval in wavelength space: lower
     
    692692     
    693693      ! All bounded intervals
    694       do jint = 2,ninterval-1
     694      DO jint = 2,ninterval-1
    695695        wavenumber1_bound = 0.01_jprb / wavelength_bound(jint)
    696696        wavenumber2_bound = 0.01_jprb / wavelength_bound(jint-1)
     
    710710      end if
    711711
    712       do jg = 1,this%ng
    713         do jin = 1,ninput
     712      DO jg = 1,this%ng
     713        DO jin = 1,ninput
    714714          mapping(jin,jg) = sum(this%gpoint_fraction(:,jg) * planck, &
    715715               &                 mask=(i_input==jin))
     
    723723
    724724      ! Loop through all intervals
    725       do jint = 1,ninterval
     725      DO jint = 1,ninterval
    726726        ! Loop through the wavenumbers for gpoint_fraction
    727         do jwav = 1,this%nwav
     727        DO jwav = 1,this%nwav
    728728          if (jint == 1) then
    729729            ! First input interval in wavelength space: lower
     
    758758      end do
    759759      if (use_fluxes_local) then
    760         do jg = 1,this%ng
     760        DO jg = 1,this%ng
    761761          mapping(:,jg) = mapping(:,jg) / sum(this%gpoint_fraction(:,jg) * planck)
    762762        end do
     
    769769    if (.not. use_fluxes_local) then
    770770      ! Normalize mapping matrix
    771       do jg = 1,size(mapping,dim=2)
     771      DO jg = 1,size(mapping,dim=2)
    772772        mapping(:,jg) = mapping(:,jg) * (1.0_jprb/sum(mapping(:,jg)))
    773773      end do
     
    828828    is_band_unassigned = .true.
    829829
    830     do jint = 1,ninterval
     830    DO jint = 1,ninterval
    831831      inext = minloc(wavelength1, dim=1, mask=is_band_unassigned)
    832832      if (jint > 1) then
     
    873873      write(nulout, '(a,i0,a,i0,a)') '  Mapping from ', nin, ' values to ', nout, ' bands (wavenumber ranges in cm-1)'
    874874      if (nout <= 40) then
    875         do jout = 1,nout
     875        DO jout = 1,nout
    876876          write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(jout)), ' to', &
    877877               &                        nint(this%wavenumber2_band(jout)), ':'
    878           do jin = 1,nin
     878          DO jin = 1,nin
    879879            write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
    880880          end do
     
    882882        end do
    883883      else
    884         do jout = 1,30
     884        DO jout = 1,30
    885885          write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(jout)), ' to', &
    886886               &                        nint(this%wavenumber2_band(jout)), ':'
    887           do jin = 1,nin
     887          DO jin = 1,nin
    888888            write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
    889889          end do
     
    893893        write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(nout)), ' to', &
    894894             &                        nint(this%wavenumber2_band(nout)), ':'
    895         do jin = 1,nin
     895        DO jin = 1,nin
    896896          write(nulout,'(f5.2)',advance='no') mapping(jin,nout)
    897897        end do
     
    901901      write(nulout, '(a,i0,a,i0,a)') '  Mapping from ', nin, ' values to ', nout, ' g-points'
    902902      if (nout <= 40) then
    903         do jout = 1,nout
     903        DO jout = 1,nout
    904904          write(nulout,'(i3,a)',advance='no') jout, ':'
    905           do jin = 1,nin
     905          DO jin = 1,nin
    906906            write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
    907907          end do
     
    909909        end do
    910910      else
    911         do jout = 1,30
     911        DO jout = 1,30
    912912          write(nulout,'(i3,a)',advance='no') jout, ':'
    913           do jin = 1,nin
     913          DO jin = 1,nin
    914914            write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
    915915          end do
     
    918918        write(nulout,'(a)') '  ...'
    919919        write(nulout,'(i3,a)',advance='no') nout, ':'
    920         do jin = 1,nin
     920        DO jin = 1,nin
    921921          write(nulout,'(f5.2)',advance='no') mapping(jin,nout)
    922922        end do
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_thermodynamics.F90

    r4773 r5158  
    143143    end if
    144144
    145     do jlev = 1,nlev
    146        do jcol = istartcol,iendcol
     145    DO jlev = 1,nlev
     146       DO jcol = istartcol,iendcol
    147147          pressure = 0.5 * (this%pressure_hl(jcol,jlev)+this%pressure_hl(jcol,jlev+1))
    148148          temperature = 0.5 * (this%temperature_hl(jcol,jlev)+this%temperature_hl(jcol,jlev+1))
     
    265265      ! be half the separation between the half-levels at the edge of
    266266      ! the two adjacent layers
    267       do jlev = 2,nlev-1
     267      DO jlev = 2,nlev-1
    268268        layer_separation(i1:i2,jlev) = (0.5_jprb * R_over_g) * temperature_hl(i1:i2,jlev+1) &
    269269             &                    * log(pressure_hl(i1:i2,jlev+2)/pressure_hl(i1:i2,jlev))
     
    276276      ! don't take the logarithm of the last pressure in each column.
    277277
    278       do jlev = 1,nlev-2
     278      DO jlev = 1,nlev-2
    279279        layer_separation(i1:i2,jlev) = (0.5_jprb * R_over_g) * temperature_hl(i1:i2,jlev+1) &
    280280             &                    * log(pressure_hl(i1:i2,jlev)/pressure_hl(i1:i2,jlev+2))
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_tripleclouds_lw.F90

    r4946 r5158  
    202202
    203203    ! Main loop over columns
    204     do jcol = istartcol, iendcol
     204    DO jcol = istartcol, iendcol
    205205      ! --------------------------------------------------------
    206206      ! Section 2: Prepare column-specific variables and arrays
     
    211211      is_clear_sky_layer = .true.
    212212      i_cloud_top = nlev+1
    213       do jlev = 1,nlev
     213      DO jlev = 1,nlev
    214214        if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
    215215          is_clear_sky_layer(jlev) = .false.
     
    268268        ! Save the spectral fluxes if required
    269269        if (config%do_save_spectral_flux) then
    270           do jlev = 1,nlev+1
     270          DO jlev = 1,nlev+1
    271271            call indexed_sum(flux_up_clear(:,jlev), &
    272272                 &           config%i_spec_from_reordered_g_lw, &
     
    294294      end if
    295295
    296       do jlev = i_cloud_top,nlev ! Start at cloud top and work down
     296      DO jlev = i_cloud_top,nlev ! Start at cloud top and work down
    297297
    298298        ! Copy over clear-sky properties
     
    308308          source_dn(:,2:,jlev)     = 0.0_jprb
    309309        else
    310           do jreg = 2,nreg
     310          DO jreg = 2,nreg
    311311            ! Cloudy sky
    312312            ! Add scaled cloud optical depth to clear-sky value
     
    358358          end do
    359359          ! Emission is scaled by the size of each region
    360           do jreg = 1,nregions
     360          DO jreg = 1,nregions
    361361            source_up(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * source_up(:,jreg,jlev)
    362362            source_dn(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * source_dn(:,jreg,jlev)
     
    375375      ! Calculate the upwelling radiation emitted by the surface, and
    376376      ! copy the surface albedo into total_albedo
    377       do jreg = 1,nregions
    378         do jg = 1,ng
     377      DO jreg = 1,nregions
     378        DO jg = 1,ng
    379379          ! region_fracs(jreg,nlev,jcol) is the fraction of each region in the
    380380          ! lowest model level
     
    387387      ! atmosphere and the total upwelling due to emission below each
    388388      ! level below using the adding method
    389       do jlev = nlev,i_cloud_top,-1
     389      DO jlev = nlev,i_cloud_top,-1
    390390
    391391        total_albedo_below        = 0.0_jprb
     
    394394          total_albedo_below = 0.0_jprb
    395395          total_source_below = 0.0_jprb
    396           do jg = 1,ng
     396          DO jg = 1,ng
    397397            inv_denom(jg,1) = 1.0_jprb &
    398398                 &  / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance(jg,1,jlev))
     
    428428          ! the operation we perform is essentially diag(total_albedo)
    429429          ! = matmul(transpose(v_matrix), diag(total_albedo_below)).
    430           do jreg = 1,nregions
    431             do jreg2 = 1,nregions
     430          DO jreg = 1,nregions
     431            DO jreg2 = 1,nregions
    432432              total_albedo(:,jreg,jlev) &
    433433                   &  = total_albedo(:,jreg,jlev) &
     
    445445      ! Section 6: Copy over downwelling fluxes above cloud top
    446446      ! --------------------------------------------------------
    447       do jlev = 1,i_cloud_top
     447      DO jlev = 1,i_cloud_top
    448448        if (config%do_clear) then
    449449          ! Clear-sky fluxes have already been averaged: use these
     
    476476             &           flux%lw_up_band(:,jcol,i_cloud_top))
    477477      end if
    478       do jlev = i_cloud_top-1,1,-1
     478      DO jlev = i_cloud_top-1,1,-1
    479479        flux_up(:,1) = trans_clear(:,jlev)*flux_up(:,1) + source_up_clear(:,jlev)
    480480        flux%lw_up(jcol,jlev) = sum(flux_up(:,1))
     
    494494      ! scattering layer, using overlap matrix to translate to the
    495495      ! regions of the first layer of cloud
    496       do jreg = 1,nregions
     496      DO jreg = 1,nregions
    497497        flux_dn(:,jreg)  = v_matrix(jreg,1,i_cloud_top,jcol)*flux_dn_clear(:,i_cloud_top)
    498498      end do
    499499
    500500      ! Final loop back down through the atmosphere to compute fluxes
    501       do jlev = i_cloud_top,nlev
     501      DO jlev = i_cloud_top,nlev
    502502
    503503        if (is_clear_sky_layer(jlev)) then
    504           do jg = 1,ng
     504          DO jg = 1,ng
    505505            flux_dn(jg,1) = (transmittance(jg,1,jlev)*flux_dn(jg,1) &
    506506                 &  + reflectance(jg,1,jlev)*total_source(jg,1,jlev+1) + source_dn(jg,1,jlev) ) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_tripleclouds_sw.F90

    r4946 r5158  
    127127    ! entire clear-sky atmosphere at once
    128128    real(jprb), dimension(config%n_g_sw, nlev) &
    129         &  :: reflectance_clear, transmittance_clear, &
     129        &  :: reflectance_clear, transmittance_clear, &
    130130         &     ref_dir_clear, trans_dir_diff_clear, trans_dir_dir_clear
    131131
     
    198198
    199199    ! Main loop over columns
    200     do jcol = istartcol, iendcol
     200    DO jcol = istartcol, iendcol
    201201      ! --------------------------------------------------------
    202202      ! Section 2: Prepare column-specific variables and arrays
     
    251251      ! cloud%crop_cloud_fraction has already been called
    252252      is_clear_sky_layer = .true.
    253       do jlev = 1,nlev
     253      DO jlev = 1,nlev
    254254        if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
    255255          is_clear_sky_layer(jlev) = .false.
     
    271271
    272272      ! Cloudy layers
    273       do jlev = 1,nlev ! Start at top-of-atmosphere
     273      DO jlev = 1,nlev ! Start at top-of-atmosphere
    274274        if (.not. is_clear_sky_layer(jlev)) then
    275           do jreg = 2,nregions
    276             do jg = 1,ng
     275          DO jreg = 2,nregions
     276            DO jg = 1,ng
    277277              ! Mapping from g-point to band
    278278              iband = config%i_band_from_reordered_g_sw(jg)
     
    315315
    316316      ! Copy surface albedos in clear-sky region
    317       do jg = 1,ng
     317      DO jg = 1,ng
    318318        total_albedo(jg,1,nlev+1) = albedo_diffuse(jg,jcol)
    319319        total_albedo_direct(jg,1,nlev+1) &
     
    324324      ! underneath
    325325      if (.not. is_clear_sky_layer(nlev)) then
    326         do jreg = 2,nregions
     326        DO jreg = 2,nregions
    327327          total_albedo(:,jreg,nlev+1)        = total_albedo(:,1,nlev+1)
    328328          total_albedo_direct(:,jreg,nlev+1) = total_albedo_direct(:,1,nlev+1)
     
    337337      ! Work up from the surface computing the total albedo of the
    338338      ! atmosphere below that point using the adding method
    339       do jlev = nlev,1,-1
     339      DO jlev = nlev,1,-1
    340340
    341341        total_albedo_below        = 0.0_jprb
     
    346346          ! "below" quantities since with no cloud overlap to worry
    347347          ! about, these are the same
    348           do jg = 1,ng
     348          DO jg = 1,ng
    349349            inv_denom(jg,1) = 1.0_jprb &
    350350                 &  / (1.0_jprb - total_albedo_clear(jg,jlev+1)*reflectance_clear(jg,jlev))
     
    361361
    362362        ! All-sky fluxes: first the clear region
    363         do jg = 1,ng
     363        DO jg = 1,ng
    364364          inv_denom(jg,1) = 1.0_jprb &
    365365               &  / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance_clear(jg,jlev))
     
    375375        ! Then the cloudy regions if any cloud is present in this layer
    376376        if (.not. is_clear_sky_layer(jlev)) then
    377           do jreg = 2,nregions
    378             do jg = 1,ng
     377          DO jreg = 2,nregions
     378            DO jg = 1,ng
    379379              inv_denom(jg,jreg) = 1.0_jprb / (1.0_jprb &
    380380                   &  - total_albedo(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev))
     
    400400          ! the operation we perform is essentially diag(total_albedo)
    401401          ! = matmul(transpose(v_matrix)), diag(total_albedo_below)).
    402           do jreg = 1,nregions
    403             do jreg2 = 1,nregions
     402          DO jreg = 1,nregions
     403            DO jreg2 = 1,nregions
    404404              total_albedo(:,jreg,jlev) &
    405405                   &  = total_albedo(:,jreg,jlev) &
     
    426426      ! Direct downwelling flux (into a plane perpendicular to the
    427427      ! sun) entering the top of each region in the top-most layer
    428       do jreg = 1,nregions
     428      DO jreg = 1,nregions
    429429        direct_dn(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol)
    430430      end do
     
    495495
    496496      ! Final loop back down through the atmosphere to compute fluxes
    497       do jlev = 1,nlev
     497      DO jlev = 1,nlev
    498498        if (config%do_clear) then
    499           do jg = 1,ng
     499          DO jg = 1,ng
    500500            flux_dn_clear(jg) = (transmittance_clear(jg,jlev)*flux_dn_clear(jg) + direct_dn_clear(jg) &
    501501               &  * (trans_dir_dir_clear(jg,jlev)*total_albedo_clear_direct(jg,jlev+1)*reflectance_clear(jg,jlev) &
     
    509509
    510510        ! All-sky fluxes: first the clear region...
    511         do jg = 1,ng
     511        DO jg = 1,ng
    512512          flux_dn(jg,1) = (transmittance_clear(jg,jlev)*flux_dn(jg,1) + direct_dn(jg,1) &
    513513               &  * (trans_dir_dir_clear(jg,jlev)*total_albedo_direct(jg,1,jlev+1)*reflectance_clear(jg,jlev) &
     
    525525          direct_dn(:,2:nregions)= 0.0_jprb
    526526        else
    527           do jreg = 2,nregions
    528             do jg = 1,ng
     527          DO jreg = 2,nregions
     528            DO jg = 1,ng
    529529              flux_dn(jg,jreg) = (transmittance(jg,jreg,jlev)*flux_dn(jg,jreg) + direct_dn(jg,jreg) &
    530530                   &  * (trans_dir_dir(jg,jreg,jlev)*total_albedo_direct(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev) &
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_two_stream.F90

    r4853 r5158  
    7272! Added for DWD (2020)
    7373!NEC$ shortloop
    74     do jg = 1, ng
     74    DO jg = 1, ng
    7575      ! Fu et al. (1997), Eq 2.9 and 2.10:
    7676      !      gamma1(jg) = LwDiffusivity * (1.0_jprb - 0.5_jprb*ssa(jg) &
     
    122122! Added for DWD (2020)
    123123!NEC$ shortloop
    124     do jg = 1, ng
     124    DO jg = 1, ng
    125125      !      gamma1(jg) = 2.0_jprb  - ssa(jg) * (1.25_jprb + 0.75_jprb*g(jg))
    126126      !      gamma2(jg) = 0.75_jprb *(ssa(jg) * (1.0_jprb - g(jg)))
     
    195195! Added for DWD (2020)
    196196!NEC$ shortloop
    197     do jg = 1, ng
     197    DO jg = 1, ng
    198198      if (od(jg) > 1.0e-3_jprd) then
    199199        k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), &
     
    291291#endif
    292292
    293     do jg = 1, ng
     293    DO jg = 1, ng
    294294      factor = (LwDiffusivityWP * 0.5_jprb) * ssa(jg)
    295295      gamma1 = LwDiffusivityWP - factor*(1.0_jprb + asymmetry(jg))
     
    379379#endif
    380380
    381     do jg = 1, ng
     381    DO jg = 1, ng
    382382      ! Compute upward and downward emission assuming the Planck
    383383      ! function to vary linearly with optical depth within the layer
     
    474474! Added for DWD (2020)
    475475!NEC$ shortloop
    476     do jg = 1, ng
     476    DO jg = 1, ng
    477477
    478478        gamma4 = 1.0_jprd - gamma3(jg)
     
    624624    trans_dir_dir = exp(trans_dir_dir)
    625625
    626     do jg = 1, ng
     626    DO jg = 1, ng
    627627
    628628      ! Zdunkowski "PIFM" (Zdunkowski et al., 1980; Contributions to
     
    652652    exponential = exp(-k_exponent*od)
    653653
    654     do jg = 1, ng
     654    DO jg = 1, ng
    655655      k_mu0 = k_exponent(jg)*mu0
    656656      one_minus_kmu0_sqr = 1.0_jprb - k_mu0*k_mu0
     
    698698#else
    699699    ! GPU-capable and vector-optimized version for ICON
    700     do jg = 1, ng
     700    DO jg = 1, ng
    701701
    702702      trans_dir_dir(jg) = max(-max(od(jg) * (1.0_jprb/mu0),0.0_jprb),-1000.0_jprb)
     
    810810! Added for DWD (2020)
    811811!NEC$ shortloop
    812     do jg = 1, ng
     812    DO jg = 1, ng
    813813      ! Note that if the minimum value is reduced (e.g. to 1.0e-24)
    814814      ! then noise starts to appear as a function of solar zenith
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/easy_netcdf.F90

    r4773 r5158  
    392392
    393393    ndimlens(:) = 0
    394     do j = 1,ndims
     394    DO j = 1,ndims
    395395      istatus = nf90_inquire_dimension(this%ncid, idimids(j), len=ndimlens(j))
    396396      if (istatus /= NF90_NOERR) then
     
    404404    if (present(ntotal)) then
    405405      ntotal = 1
    406       do j = 1, ndims
     406      DO j = 1, ndims
    407407        ntotal = ntotal * ndimlens(j)
    408408      end do
     
    613613    ! Compute number of elements of the variable in the file
    614614    ntotal = 1
    615     do j = 1, ndims
     615    DO j = 1, ndims
    616616      ntotal = ntotal * ndimlens(j)
    617617    end do
     
    658658    ! Compute number of elements of the variable in the file
    659659    ntotal = 1
    660     do j = 1, ndims
     660    DO j = 1, ndims
    661661      ntotal = ntotal * ndimlens(j)
    662662    end do
     
    707707    ! i.e. excluding the slowest varying dimension, indexed by "index"
    708708    ntotal = 1
    709     do j = 1, ndims-1
     709    DO j = 1, ndims-1
    710710      ntotal = ntotal * ndimlens(j)
    711711    end do
     
    761761    ! Ensure variable has only one dimension in the file
    762762    n = 1
    763     do j = 1, ndims
     763    DO j = 1, ndims
    764764      n = n * ndimlens(j)
    765765      if (j > 1 .and. ndimlens(j) > 1) then
     
    819819    ! Ensure variable has only one dimension in the file
    820820    n = 1
    821     do j = 1, ndims
     821    DO j = 1, ndims
    822822      n = n * ndimlens(j)
    823823      if (j > 1 .and. ndimlens(j) > 1) then
     
    878878    ! Ensure variable has only one dimension in the file
    879879    n = 1
    880     do j = 1, ndims
     880    DO j = 1, ndims
    881881      n = n * ndimlens(j)
    882882      if (j > 1 .and. ndimlens(j) > 1) then
     
    938938    ! Ensure variable has only one dimension aside from the last one
    939939    n = 1
    940     do j = 1, ndims-1
     940    DO j = 1, ndims-1
    941941      n = n * ndimlens(j)
    942942      if (j > 1 .and. ndimlens(j) > 1) then
     
    10211021    ! dimensions
    10221022    ntotal = 1
    1023     do j = 1, ndims
     1023    DO j = 1, ndims
    10241024      ntotal = ntotal * ndimlens(j)
    10251025      if (j > 2 .and. ndimlens(j) > 1) then
     
    11331133    ! dimensions
    11341134    ntotal = 1
    1135     do j = 1, ndims
     1135    DO j = 1, ndims
    11361136      ntotal = ntotal * ndimlens(j)
    11371137      if (j > 2 .and. ndimlens(j) > 1) then
     
    12021202      vcount(1:2) = [ndimlen1,1]
    12031203     
    1204       do j = 1,ndimlen2
     1204      DO j = 1,ndimlen2
    12051205        vstart(2) = j
    12061206        istatus = nf90_get_var(this%ncid, ivarid, matrix(:,j), start=vstart, count=vcount)
     
    12521252    ! dimensions aside from the last one
    12531253    ntotal = 1
    1254     do j = 1, ndims-1
     1254    DO j = 1, ndims-1
    12551255      ntotal = ntotal * ndimlens(j)
    12561256      if (j > 2 .and. ndimlens(j) > 1) then
     
    13761376    ! dimensions
    13771377    ntotal = 1
    1378     do j = 1, ndims
     1378    DO j = 1, ndims
    13791379      ntotal = ntotal * ndimlens(j)
    13801380      if (j > 3 .and. ndimlens(j) > 1) then
     
    15121512    ! dimensions aside from the last one
    15131513    ntotal = 1
    1514     do j = 1, ndims-1
     1514    DO j = 1, ndims-1
    15151515      ntotal = ntotal * ndimlens(j)
    15161516      if (j > 3 .and. ndimlens(j) > 1) then
     
    16541654    ! dimensions
    16551655    ntotal = 1
    1656     do j = 1, ndims
     1656    DO j = 1, ndims
    16571657      ntotal = ntotal * ndimlens(j)
    16581658      if (j > 4 .and. ndimlens(j) > 1) then
     
    18041804
    18051805    ! Pad with blanks since nf90_get_att does not do this
    1806     do j = i_attr_len+1,len(attr_str)
     1806    DO j = i_attr_len+1,len(attr_str)
    18071807      attr_str(j:j) = ' '
    18081808    end do
     
    18781878
    18791879    ! Pad with blanks since nf90_get_att does not do this
    1880     do j = i_attr_len+1,len(attr_str)
     1880    DO j = i_attr_len+1,len(attr_str)
    18811881      attr_str(j:j) = ' '
    18821882    end do
     
    27052705    end if
    27062706
    2707     do jdim = 1,ndims
     2707    DO jdim = 1,ndims
    27082708      istatus = nf90_inquire_dimension(infile%ncid, idimids(jdim), &
    27092709           &  name=dimname, len=dimlen)
     
    27712771
    27722772    ! Map dimension IDs
    2773     do jdim = 1,ndims
     2773    DO jdim = 1,ndims
    27742774      istatus = nf90_inquire_dimension(infile%ncid, idimids_in(jdim), name=dim_name)
    27752775      if (istatus /= NF90_NOERR) then
     
    28012801
    28022802    ! Copy attributes
    2803     do jattr = 1,nattr
     2803    DO jattr = 1,nattr
    28042804      istatus = nf90_inq_attname(infile%ncid, ivarid_in, jattr, attr_name)
    28052805      if (istatus /= NF90_NOERR) then
  • LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/print_matrix.F90

    r4773 r5158  
    1616      write(unit_local,'(a,a,$)') name, '=['
    1717    end if
    18     do i = 1,size(vec,1)
     18    DO i = 1,size(vec,1)
    1919       write(unit_local,'(f16.8,$)') vec(i)
    2020    end do
     
    4343      write(unit_local,'(a,a,$)') name, '=['
    4444    end if
    45     do i = 1,size(mat,1)
    46        do j = 1,size(mat,2)
     45    DO i = 1,size(mat,1)
     46       DO j = 1,size(mat,2)
    4747          write(unit_local,'(f16.8,$)') mat(i,j)
    4848       end do
     
    6161    integer    :: i, j, k
    6262    write(6,'(a,a)') name, '='
    63     do k = 1,size(mat,3)
    64       do i = 1,size(mat,1)
    65         do j = 1,size(mat,2)
     63    DO k = 1,size(mat,3)
     64      DO i = 1,size(mat,1)
     65        DO j = 1,size(mat,2)
    6666             write(6,'(f16.8,$)') mat(i,j,k)
    6767          end do
Note: See TracChangeset for help on using the changeset viewer.