Changeset 5158 for LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad
- Timestamp:
- Aug 2, 2024, 2:12:03 PM (6 months ago)
- 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 311 311 tstart = omp_get_wtime() 312 312 #endif 313 dojrepeat = 1,driver_config%nrepeat313 DO jrepeat = 1,driver_config%nrepeat 314 314 315 315 if (driver_config%do_parallel) then … … 321 321 322 322 !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME) 323 dojblock = 1, nblock323 DO jblock = 1, nblock 324 324 ! Specify the range of columns to process. 325 325 istartcol = (jblock-1) * driver_config%nblocksize & -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_config.F90
r4773 r5158 400 400 ! wavelength bound, noting that the number of diagnostics is one 401 401 ! fewer than the number of valid bounds 402 dojdiag = 0,NMaxSpectralDiag402 DO jdiag = 0,NMaxSpectralDiag 403 403 if (this%sw_diag_wavelength_bound(jdiag+1) < 0.0_jprb) then 404 404 this%n_sw_diag = max(0,jdiag-1) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_driver_read_input.F90
r4773 r5158 563 563 564 564 ! Loop through all radiatively important gases 565 dojgas = 1,NMaxGases565 DO jgas = 1,NMaxGases 566 566 if (jgas == IH2O) then 567 567 if (file%exists('q')) then -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver.F90
r4773 r5158 369 369 tstart = omp_get_wtime() 370 370 #endif 371 dojrepeat = 1,driver_config%nrepeat371 DO jrepeat = 1,driver_config%nrepeat 372 372 373 373 ! if (driver_config%do_parallel) then … … 379 379 380 380 !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME) 381 dojblock = 1, nblock381 DO jblock = 1, nblock 382 382 ! Specify the range of columns to process. 383 383 istartcol = (jblock-1) * driver_config%nblocksize & -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ecrad_ifs_driver_blocked.F90
r4773 r5158 371 371 tstart = omp_get_wtime() 372 372 #endif 373 dojrepeat = 1,driver_config%nrepeat373 DO jrepeat = 1,driver_config%nrepeat 374 374 375 375 ! if (driver_config%do_parallel) then … … 378 378 !$OMP PARALLEL DO SCHEDULE(DYNAMIC,1)& 379 379 !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB) 380 dojrl=1,ncol,nproma380 DO jrl=1,ncol,nproma 381 381 ibeg=jrl 382 382 iend=min(ibeg+nproma-1,ncol) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/ifs_blocking.F90
r4773 r5158 329 329 !$OMP PARALLEL DO SCHEDULE(RUNTIME)& 330 330 !$OMP&PRIVATE(IB,IFLD) 331 doib=1,ngpblks332 doifld=1,ifs_config%ifldstot331 DO ib=1,ngpblks 332 DO ifld=1,ifs_config%ifldstot 333 333 zrgp(:,ifld,ib) = 0._jprb 334 334 enddo … … 351 351 !$OMP PARALLEL DO SCHEDULE(RUNTIME)& 352 352 !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB,JAER,JOFF,JLEV,JALB) 353 dojrl=1,ncol,nproma353 DO jrl=1,ncol,nproma 354 354 355 355 ibeg=jrl … … 363 363 zrgp(1:il,ifs_config%iamu0,ib) = single_level%cos_sza(ibeg:iend) ! cosine of solar zenith ang (mu0) 364 364 365 dojemiss=1,yderad%nlwemiss365 DO jemiss=1,yderad%nlwemiss 366 366 zrgp(1:il,ifs_config%iemiss+jemiss-1,ib) = single_level%lw_emissivity(ibeg:iend,jemiss) 367 367 enddo … … 378 378 ! zrgp(1:il,islon,ib) = ??? 379 379 380 dojalb=1,yderad%nsw380 DO jalb=1,yderad%nsw 381 381 zrgp(1:il,ifs_config%iald+jalb-1,ib) = single_level%sw_albedo(ibeg:iend,jalb) 382 382 enddo 383 383 384 384 if (allocated(single_level%sw_albedo_direct)) then 385 dojalb=1,yderad%nsw385 DO jalb=1,yderad%nsw 386 386 zrgp(1:il,ifs_config%ialp+jalb-1,ib) = single_level%sw_albedo_direct(ibeg:iend,jalb) 387 387 end do 388 388 else 389 dojalb=1,yderad%nsw389 DO jalb=1,yderad%nsw 390 390 zrgp(1:il,ifs_config%ialp+jalb-1,ib) = single_level%sw_albedo(ibeg:iend,jalb) 391 391 end do 392 392 end if 393 393 394 dojlev=1,nlev394 DO jlev=1,nlev 395 395 zrgp(1:il,ifs_config%iti+jlev-1,ib) = temperature_fl(ibeg:iend,jlev) ! full level temperature 396 396 zrgp(1:il,ifs_config%ipr+jlev-1,ib) = pressure_fl(ibeg:iend,jlev) ! full level pressure … … 398 398 enddo 399 399 400 dojlev=1,nlev400 DO jlev=1,nlev 401 401 zrgp(1:il,ifs_config%iwv+jlev-1,ib) = gas%mixing_ratio(ibeg:iend,jlev,IH2O) ! this is already in MassMixingRatio units 402 402 if (rad_config%do_clouds) then … … 421 421 if (yderad%naermacc == 1) then 422 422 joff=ifs_config%iaero 423 dojaer=1,rad_config%n_aerosol_types424 dojlev=1,nlev423 DO jaer=1,rad_config%n_aerosol_types 424 DO jlev=1,nlev 425 425 zrgp(1:il,joff,ib) = aerosol%mixing_ratio(ibeg:iend,jlev,jaer) 426 426 joff=joff+1 … … 429 429 endif 430 430 431 dojlev=1,nlev+1431 DO jlev=1,nlev+1 432 432 ! zrgp(1:il,ihpr+jlev-1,ib) = ??? 433 433 zrgp(1:il,ifs_config%iaprs+jlev-1,ib) = thermodynamics%pressure_hl(ibeg:iend,jlev) … … 453 453 ! local workaround variables for standalone input files 454 454 if (rad_config%do_clouds) then 455 dojlev=1,nlev455 DO jlev=1,nlev 456 456 ! missing full-level temperature and pressure as well as land-sea-mask 457 457 zrgp(1:il,ifs_config%ire_liq+jlev-1,ib) = cloud%re_liq(ibeg:iend,jlev) 458 458 zrgp(1:il,ifs_config%ire_ice+jlev-1,ib) = cloud%re_ice(ibeg:iend,jlev) 459 459 enddo 460 dojlev=1,nlev-1460 DO jlev=1,nlev-1 461 461 ! for the love of it, I can't figure this one out. Probably to do with 462 462 ! my crude approach of setting PGEMU? … … 465 465 if(present(iseed)) iseed(1:il,ib) = single_level%iseed(ibeg:iend) 466 466 else 467 dojlev=1,nlev467 DO jlev=1,nlev 468 468 ! missing full-level temperature and pressure as well as land-sea-mask 469 469 zrgp(1:il,ifs_config%ire_liq+jlev-1,ib) = 0._jprb 470 470 zrgp(1:il,ifs_config%ire_ice+jlev-1,ib) = 0._jprb 471 471 enddo 472 dojlev=1,nlev-1472 DO jlev=1,nlev-1 473 473 zrgp(1:il,ifs_config%ioverlap+jlev-1,ib) = 0._jprb 474 474 enddo … … 531 531 !$OMP PARALLEL DO SCHEDULE(RUNTIME)& 532 532 !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB,JLEV,JG) 533 dojrl=1,ncol,nproma533 DO jrl=1,ncol,nproma 534 534 ibeg=jrl 535 535 iend=min(ibeg+nproma-1,ncol) … … 537 537 ib=(jrl-1)/nproma+1 538 538 539 dojlev=1,nlev+1539 DO jlev=1,nlev+1 540 540 flux%sw_up(ibeg:iend,jlev) = zrgp(1:il,ifs_config%ifrso+jlev-1,ib) 541 541 flux%lw_up(ibeg:iend,jlev) = zrgp(1:il,ifs_config%ifrth+jlev-1,ib) … … 561 561 emissivity_out(ibeg:iend) = zrgp(1:il,ifs_config%iemit,ib) 562 562 if (yradiation%yrerad%lapproxswupdate) then 563 dojg=1,yradiation%yrerad%nsw563 DO jg=1,yradiation%yrerad%nsw 564 564 flux_diffuse_band(ibeg:iend,jg) = zrgp(1:il,ifs_config%iswdiffuseband+jg-1,ib) 565 565 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 57 57 & od_scaling, total_cloud_cover) 58 58 59 dojlev = 1,nlev60 dojcol = 1,ncol59 DO jlev = 1,nlev 60 DO jcol = 1,ncol 61 61 write(6,'(f5.2,a)',advance='no') od_scaling(jcol,jlev), ' ' 62 62 end do -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_random_number_generator.F90
r4773 r5158 29 29 call random_number_generator%initialize(itype=IRngType, iseed=212075152, & 30 30 & nmaxstreams=streammax) 31 dojl = 1,131 DO jl = 1,1 32 32 print *, 'initial_state = [ ', int(random_number_generator%istate(1:streammax),jpib), ' ]' 33 33 call random_number_generator%uniform_distribution(vec) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/driver/test_spartacus_math.F90
r4773 r5158 45 45 46 46 write(*,*) 'A =' 47 doj = 1,m47 DO j = 1,m 48 48 write(*,*) A(1,j,:) 49 49 end do 50 50 51 51 write(*,*) 'B =' 52 doj = 1,m52 DO j = 1,m 53 53 write(*,*) B(1,j,:) 54 54 end do … … 60 60 61 61 write(*,*) 'C = A*B =' 62 doj = 1,m62 DO j = 1,m 63 63 write(*,*) C(1,j,:) 64 64 end do … … 70 70 71 71 write(*,*) 'C = A\B =' 72 doj = 1,m72 DO j = 1,m 73 73 write(*,*) C(1,j,:) 74 74 end do … … 78 78 79 79 write(*,*) 'expm(A) =' 80 doj = 1,m80 DO j = 1,m 81 81 write(*,*) A(1,j,:) 82 82 end do … … 93 93 94 94 write(*,*) 'fast_expm(A) = ' 95 doj = 1,m95 DO j = 1,m 96 96 write(*,*) A(1,j,:) 97 97 end do … … 106 106 107 107 write(*,*) 'expm(A) = ' 108 doj = 1,m108 DO j = 1,m 109 109 write(*,*) A(1,j,:) 110 110 end do … … 115 115 116 116 write(*,*) 'expm(zeros) = ' 117 doj = 1,m117 DO j = 1,m 118 118 write(*,*) A(1,j,:) 119 119 end do … … 124 124 125 125 write(*,*) 'fast_expm(A) = ' 126 doj = 1,m126 DO j = 1,m 127 127 write(*,*) A(1,j,:) 128 128 end do … … 140 140 141 141 write(*,*) 'expm(A) = ' 142 doj = 1,m142 DO j = 1,m 143 143 write(*,*) A(1,j,:) 144 144 end do … … 149 149 150 150 write(*,*) 'expm(zeros) = ' 151 doj = 1,m151 DO j = 1,m 152 152 write(*,*) A(1,j,:) 153 153 end do -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/ifs/fcttre.func.h
r4876 r5158 170 170 FOEELIQ2ICE( PTARE ) = EXP( (PTARE-RTT)*(R3LES/(PTARE-R4LES) - R3IES/(PTARE-R4IES))) 171 171 FOELSON( PTARE ) = EXP( -6096.9385_JPRB/PTARE + 21.2409642_JPRB & 172 172 - 2.711193E-2_JPRB * PTARE & 173 173 + 1.673952E-5_JPRB * PTARE**2 & 174 174 + 2.433502_JPRB * LOG(PTARE)) 175 175 176 176 REAL(KIND=JPRB) :: FOEEWM_V,FOEEWMCU_V,FOELES_V,FOEIES_V -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/lmdz/aeropt_5wv_ecrad.F90
r4853 r5158 54 54 55 55 ! Loop over level 56 dojlev = istartlev,iendlev56 DO jlev = istartlev,iendlev 57 57 58 58 ! Loop over column 59 dojcol = istartcol,iendcol59 DO jcol = istartcol,iendcol 60 60 61 61 ! Compute relative humidity with respect to liquid … … 71 71 od_aerosol_mono = 0.0_jprb 72 72 73 dojtype = 1,config%n_aerosol_types73 DO jtype = 1,config%n_aerosol_types 74 74 ! Add the optical depth for this aerosol type to the total for all aerosols. 75 75 ! Note that the following expressions are array-wise -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_lw.F90
r4773 r5158 83 83 ! also the "source", which is the upwelling flux due to emission 84 84 ! below that level 85 dojlev = nlev,1,-185 DO jlev = nlev,1,-1 86 86 ! Next loop over columns. We could do this by indexing the 87 87 ! entire inner dimension as follows, e.g. for the first line: … … 90 90 ! routine by a factor of 2! Rather, we do it with an explicit 91 91 ! loop. 92 dojcol = 1,ncol92 DO jcol = 1,ncol 93 93 ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10: 94 94 inv_denominator(jcol,jlev) = 1.0_jprb & … … 114 114 ! Work back down through the atmosphere computing the fluxes at 115 115 ! each half-level 116 dojlev = 1,nlev117 dojcol = 1,ncol116 DO jlev = 1,nlev 117 DO jcol = 1,ncol 118 118 ! Shonk & Hogan (2008) Eq 14 (after simplification): 119 119 flux_dn(jcol,jlev+1) & … … 201 201 ! also the "source", which is the upwelling flux due to emission 202 202 ! below that level 203 dojlev = nlev,i_cloud_top,-1203 DO jlev = nlev,i_cloud_top,-1 204 204 if (is_clear_sky_layer(jlev)) then 205 205 ! Reflectance of this layer is zero, simplifying the expression 206 dojcol = 1,ncol206 DO jcol = 1,ncol 207 207 albedo(jcol,jlev) = transmittance(jcol,jlev)*transmittance(jcol,jlev)*albedo(jcol,jlev+1) 208 208 source(jcol,jlev) = source_up(jcol,jlev) & … … 212 212 else 213 213 ! Loop over columns; explicit loop seems to be faster 214 dojcol = 1,ncol214 DO jcol = 1,ncol 215 215 ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10: 216 216 inv_denominator(jcol,jlev) = 1.0_jprb & … … 231 231 flux_up(:,i_cloud_top) = source(:,i_cloud_top) & 232 232 & + albedo(:,i_cloud_top)*flux_dn(:,i_cloud_top) 233 dojlev = i_cloud_top-1,1,-1233 DO jlev = i_cloud_top-1,1,-1 234 234 flux_up(:,jlev) = transmittance(:,jlev)*flux_up(:,jlev+1) + source_up(:,jlev) 235 235 end do … … 237 237 ! Work back down through the atmosphere from cloud top computing 238 238 ! the fluxes at each half-level 239 dojlev = i_cloud_top,nlev239 DO jlev = i_cloud_top,nlev 240 240 if (is_clear_sky_layer(jlev)) then 241 dojcol = 1,ncol241 DO jcol = 1,ncol 242 242 flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) & 243 243 & + source_dn(jcol,jlev) … … 246 246 end do 247 247 else 248 dojcol = 1,ncol248 DO jcol = 1,ncol 249 249 ! Shonk & Hogan (2008) Eq 14 (after simplification): 250 250 flux_dn(jcol,jlev+1) & … … 309 309 ! Added for DWD (2020) 310 310 !NEC$ outerloop_unroll(8) 311 dojlev = 1,nlev312 dojcol = 1,ncol311 DO jlev = 1,nlev 312 DO jcol = 1,ncol 313 313 flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) + source_dn(jcol,jlev) 314 314 end do … … 322 322 ! Added for DWD (2020) 323 323 !NEC$ outerloop_unroll(8) 324 dojlev = nlev,1,-1325 dojcol = 1,ncol324 DO jlev = nlev,1,-1 325 DO jcol = 1,ncol 326 326 flux_up(jcol,jlev) = transmittance(jcol,jlev)*flux_up(jcol,jlev+1) + source_up(jcol,jlev) 327 327 end do -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_adding_ica_sw.F90
r4773 r5158 84 84 ! half-level by working down through the atmosphere 85 85 flux_dn_direct(:,1) = incoming_toa 86 dojlev = 1,nlev86 DO jlev = 1,nlev 87 87 flux_dn_direct(:,jlev+1) = flux_dn_direct(:,jlev)*trans_dir_dir(:,jlev) 88 88 end do … … 100 100 ! Added for DWD (2020) 101 101 !NEC$ outerloop_unroll(8) 102 dojlev = nlev,1,-1102 DO jlev = nlev,1,-1 103 103 ! Next loop over columns. We could do this by indexing the 104 104 ! entire inner dimension as follows, e.g. for the first line: … … 107 107 ! routine by a factor of 2! Rather, we do it with an explicit 108 108 ! loop. 109 dojcol = 1,ncol109 DO jcol = 1,ncol 110 110 ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10: 111 111 inv_denominator(jcol,jlev) = 1.0_jprb / (1.0_jprb-albedo(jcol,jlev+1)*reflectance(jcol,jlev)) … … 132 132 ! Added for DWD (2020) 133 133 !NEC$ outerloop_unroll(8) 134 dojlev = 1,nlev135 dojcol = 1,ncol134 DO jlev = 1,nlev 135 DO jcol = 1,ncol 136 136 ! Shonk & Hogan (2008) Eq 14 (after simplification): 137 137 flux_dn_diffuse(jcol,jlev+1) & -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics.F90
r4853 r5158 243 243 244 244 if (ao%use_hydrophilic) then 245 dojtype = 1,n_type_philic245 DO jtype = 1,n_type_philic 246 246 ao%mass_ext_sw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype)) 247 247 ao%ssa_sw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype) & … … 266 266 267 267 if (ao%use_hydrophilic) then 268 dojtype = 1,n_type_philic268 DO jtype = 1,n_type_philic 269 269 ao%mass_ext_lw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype)) 270 270 ao%ssa_lw_philic(:,:,jtype) = matmul(mapping, mass_ext_philic(:,:,jtype) & … … 280 280 if (allocated(ao%wavelength_mono)) then 281 281 ! Monochromatic wavelengths also required 282 dojwl = 1,nmono282 DO jwl = 1,nmono 283 283 ! Wavelength (m) to wavenumber (cm-1) 284 284 wavenumber_target = 0.01_jprb / ao%wavelength_mono(jwl) … … 292 292 else 293 293 iwn = 1 294 dowhile (wavenumber(iwn+1) < wavenumber_target .and. iwn < nwn-1)294 DO while (wavenumber(iwn+1) < wavenumber_target .and. iwn < nwn-1) 295 295 iwn = iwn + 1 296 296 end do … … 419 419 420 420 if (ao%use_hydrophilic) then 421 dojtype = 1,ao%n_type_philic421 DO jtype = 1,ao%n_type_philic 422 422 ao%mass_ext_sw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_sw_philic(:,:,jtype)) 423 423 ao%ssa_sw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_sw_philic(:,:,jtype) & … … 451 451 452 452 if (ao%use_hydrophilic) then 453 dojtype = 1,ao%n_type_philic453 DO jtype = 1,ao%n_type_philic 454 454 ao%mass_ext_lw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_lw_philic(:,:,jtype)) 455 455 ao%ssa_lw_philic(:,:,jtype) = matmul(mapping, ao_legacy%mass_ext_lw_philic(:,:,jtype) & … … 579 579 ! Aerosol mixing ratios have been provided 580 580 581 dojtype = 1,config%n_aerosol_types581 DO jtype = 1,config%n_aerosol_types 582 582 if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then 583 583 write(nulerr,'(a)') '*** Error: not all aerosol types are defined' … … 612 612 613 613 ! Loop over column 614 dojcol = istartcol,iendcol614 DO jcol = istartcol,iendcol 615 615 616 616 ! Reset temporary arrays … … 622 622 scat_g_lw_aerosol = 0.0_jprb 623 623 624 dojlev = istartlev,iendlev624 DO jlev = istartlev,iendlev 625 625 ! Compute relative humidity with respect to liquid 626 626 ! saturation and the index to the relative-humidity index of … … 634 634 end do 635 635 636 dojtype = 1,config%n_aerosol_types636 DO jtype = 1,config%n_aerosol_types 637 637 itype = ao%itype(jtype) 638 638 … … 643 643 ! dimension being spectral band. 644 644 if (ao%iclass(jtype) == IAerosolClassHydrophobic) then 645 dojlev = istartlev,iendlev645 DO jlev = istartlev,iendlev 646 646 mixing_ratio = aerosol%mixing_ratio(jcol,jlev,jtype) 647 dojband = 1,config%n_bands_sw647 DO jband = 1,config%n_bands_sw 648 648 local_od_sw = factor(jlev) * mixing_ratio & 649 649 & * ao%mass_ext_sw_phobic(jband,itype) … … 656 656 end do 657 657 if (config%do_lw_aerosol_scattering) then 658 dojband = 1,config%n_bands_lw658 DO jband = 1,config%n_bands_lw 659 659 local_od_lw = factor(jlev) * mixing_ratio & 660 660 & * ao%mass_ext_lw_phobic(jband,itype) … … 670 670 ! weight the optical depth by the single scattering 671 671 ! co-albedo 672 dojband = 1,config%n_bands_lw672 DO jband = 1,config%n_bands_lw 673 673 od_lw_aerosol(jband,jlev) = od_lw_aerosol(jband,jlev) & 674 674 & + factor(jlev) * mixing_ratio & … … 682 682 ! Hydrophilic aerosols require the look-up tables to 683 683 ! be indexed with irh 684 dojlev = istartlev,iendlev684 DO jlev = istartlev,iendlev 685 685 mixing_ratio = aerosol%mixing_ratio(jcol,jlev,jtype) 686 686 irh = irhs(jlev) 687 dojband = 1,config%n_bands_sw687 DO jband = 1,config%n_bands_sw 688 688 local_od_sw = factor(jlev) * mixing_ratio & 689 689 & * ao%mass_ext_sw_philic(jband,irh,itype) … … 696 696 end do 697 697 if (config%do_lw_aerosol_scattering) then 698 dojband = 1,config%n_bands_lw698 DO jband = 1,config%n_bands_lw 699 699 local_od_lw = factor(jlev) * mixing_ratio & 700 700 & * ao%mass_ext_lw_philic(jband,irh,itype) … … 710 710 ! weight the optical depth by the single scattering 711 711 ! co-albedo 712 dojband = 1,config%n_bands_lw712 DO jband = 1,config%n_bands_lw 713 713 od_lw_aerosol(jband,jlev) = od_lw_aerosol(jband,jlev) & 714 714 & + factor(jlev) * mixing_ratio & … … 740 740 741 741 ! We can assume the band and g-point indices are the same 742 dojlev = 1,nlev743 dojg = 1,config%n_g_sw742 DO jlev = 1,nlev 743 DO jg = 1,config%n_g_sw 744 744 local_scat = ssa_sw(jg,jlev,jcol)*od_sw(jg,jlev,jcol) + scat_sw_aerosol(jg,jlev) 745 745 od_sw(jg,jlev,jcol) = od_sw(jg,jlev,jcol) + od_sw_aerosol(jg,jlev) … … 751 751 else 752 752 753 dojlev = 1,nlev754 dojg = 1,config%n_g_sw753 DO jlev = 1,nlev 754 DO jg = 1,config%n_g_sw 755 755 ! Need to map between bands and g-points 756 756 iband = config%i_band_from_reordered_g_sw(jg) … … 781 781 & scat_lw_aerosol, scat_g_lw_aerosol) 782 782 783 dojlev = istartlev,iendlev784 dojg = 1,config%n_g_lw783 DO jlev = istartlev,iendlev 784 DO jg = 1,config%n_g_lw 785 785 iband = config%i_band_from_reordered_g_lw(jg) 786 786 local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband,jlev) … … 802 802 if (config%do_cloud_aerosol_per_lw_g_point) then 803 803 ! We can assume band and g-point indices are the same 804 dojlev = istartlev,iendlev805 dojg = 1,config%n_g_lw804 DO jlev = istartlev,iendlev 805 DO jg = 1,config%n_g_lw 806 806 od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) + od_lw_aerosol(jg,jlev) 807 807 end do 808 808 end do 809 809 else 810 dojlev = istartlev,iendlev811 dojg = 1,config%n_g_lw810 DO jlev = istartlev,iendlev 811 DO jg = 1,config%n_g_lw 812 812 od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) & 813 813 & + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg),jlev) … … 901 901 902 902 ! Loop over position 903 dojcol = istartcol,iendcol903 DO jcol = istartcol,iendcol 904 904 ! Added for DWD (2020) 905 905 !NEC$ forced_collapse 906 dojlev = istartlev,iendlev907 dojb = 1,config%n_bands_sw906 DO jlev = istartlev,iendlev 907 DO jb = 1,config%n_bands_sw 908 908 od_sw_aerosol(jb,jlev) = aerosol%od_sw(jb,jlev,jcol) 909 909 scat_sw_aerosol(jb,jlev) = aerosol%ssa_sw(jb,jlev,jcol) * od_sw_aerosol(jb,jlev) … … 923 923 ! properties (noting that any gas scattering will have an 924 924 ! asymmetry factor of zero) 925 dojlev = istartlev,iendlev925 DO jlev = istartlev,iendlev 926 926 if (od_sw_aerosol(1,jlev) > 0.0_jprb) then 927 dojg = 1,config%n_g_sw927 DO jg = 1,config%n_g_sw 928 928 iband = config%i_band_from_reordered_g_sw(jg) 929 929 local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband,jlev) … … 961 961 962 962 ! Loop over position 963 dojcol = istartcol,iendcol963 DO jcol = istartcol,iendcol 964 964 ! Added for DWD (2020) 965 965 !NEC$ forced_collapse 966 dojlev = istartlev,iendlev967 dojb = 1,config%n_bands_lw966 DO jlev = istartlev,iendlev 967 DO jb = 1,config%n_bands_lw 968 968 od_lw_aerosol(jb,jlev) = aerosol%od_lw(jb,jlev,jcol) 969 969 scat_lw_aerosol(jb,jlev) = aerosol%ssa_lw(jb,jlev,jcol) * od_lw_aerosol(jb,jlev) … … 974 974 end do 975 975 end do 976 dojlev = istartlev,iendlev977 dojg = 1,config%n_g_lw976 DO jlev = istartlev,iendlev 977 DO jg = 1,config%n_g_lw 978 978 iband = config%i_band_from_reordered_g_lw(jg) 979 979 if (od_lw_aerosol(iband,jlev) > 0.0_jprb) then … … 995 995 996 996 ! Loop over position 997 dojcol = istartcol,iendcol997 DO jcol = istartcol,iendcol 998 998 ! Added for DWD (2020) 999 999 !NEC$ forced_collapse 1000 dojlev = istartlev,iendlev1000 DO jlev = istartlev,iendlev 1001 1001 ! If aerosol longwave scattering is not included then we 1002 1002 ! weight the optical depth by the single scattering 1003 1003 ! co-albedo 1004 dojb = 1, config%n_bands_lw1004 DO jb = 1, config%n_bands_lw 1005 1005 od_lw_aerosol(jb,jlev) = aerosol%od_lw(jb,jlev,jcol) & 1006 1006 & * (1.0_jprb - aerosol%ssa_lw(jb,jlev,jcol)) 1007 1007 end do 1008 1008 end do 1009 dojlev = istartlev,iendlev1010 dojg = 1,config%n_g_lw1009 DO jlev = istartlev,iendlev 1010 DO jg = 1,config%n_g_lw 1011 1011 od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) & 1012 1012 & + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg),jlev) … … 1120 1120 if (lhook) call dr_hook('radiation_aerosol_optics:aerosol_extinction',0,hook_handle) 1121 1121 1122 dojtype = 1,config%n_aerosol_types1122 DO jtype = 1,config%n_aerosol_types 1123 1123 if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then 1124 1124 write(nulerr,'(a)') '*** Error: not all aerosol types are defined' … … 1138 1138 1139 1139 ! Loop over position 1140 dojcol = istartcol,iendcol1140 DO jcol = istartcol,iendcol 1141 1141 ext = 0.0_jprb 1142 1142 ! Get relative-humidity index 1143 1143 irh = ao%calc_rh_index(relative_humidity(jcol)) 1144 1144 ! Add extinction coefficients from each aerosol type 1145 dojtype = 1,config%n_aerosol_types1145 DO jtype = 1,config%n_aerosol_types 1146 1146 if (ao%iclass(jtype) == IAerosolClassHydrophobic) then 1147 1147 ext = ext + mixing_ratio(jcol,jtype) & -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_data.F90
r4860 r5158 619 619 iend = ubound(itypes,1) 620 620 621 dojtype = istart, iend621 DO jtype = istart, iend 622 622 if (itypes(jtype) == 0) then 623 623 call this%set_empty_type(jtype) … … 655 655 else 656 656 calc_rh_index = 1 657 dowhile (rh > this%rh_lower(calc_rh_index + 1))657 DO while (rh > this%rh_lower(calc_rh_index + 1)) 658 658 calc_rh_index = calc_rh_index + 1 659 659 end do … … 682 682 end if 683 683 684 dojtype = 1,size(i_type_map)684 DO jtype = 1,size(i_type_map) 685 685 if (i_type_map(jtype) > 0) then 686 686 write(nulout,'(i4,a,a)') jtype, ' -> hydrophobic type ', & … … 714 714 715 715 ! Find index of first character 716 dowhile (i_line_current < iline)716 DO while (i_line_current < iline) 717 717 i_start_new = scan(str(istart:iend), new_line(' ')) 718 718 if (i_start_new == 0) then -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_description.F90
r4853 r5158 165 165 166 166 ! Loop over hydrophilic types 167 doja = 1,size(this%bin_philic)167 DO ja = 1,size(this%bin_philic) 168 168 ! Check if we have a match 169 169 if (to_string(this%code_philic(:,ja)) == code_str & … … 176 176 end do 177 177 ! Repeat for the hydrophobic types 178 doja = 1,size(this%bin_phobic)178 DO ja = 1,size(this%bin_phobic) 179 179 if (to_string(this%code_phobic(:,ja)) == code_str & 180 180 & .and. trim(to_string(this%optical_model_phobic(:,ja))) & … … 255 255 if (lhydrophilic) then 256 256 ! Loop over hydrophilic aerosol types 257 doja = 1,size(this%bin_philic)257 DO ja = 1,size(this%bin_philic) 258 258 current_score = 0 259 259 if (to_string(this%code_philic(:,ja)) == code_str) then … … 306 306 else 307 307 ! Loop over hydrophobic aerosol types 308 doja = 1,size(this%bin_phobic)308 DO ja = 1,size(this%bin_phobic) 309 309 current_score = 0 310 310 if (to_string(this%code_phobic(:,ja)) == code_str) then … … 376 376 character(len=size(arr)) :: str 377 377 integer :: jc 378 dojc = 1,size(arr)378 DO jc = 1,size(arr) 379 379 if (ichar(arr(jc)) == 0) then 380 380 ! Replace NULL character with a space -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud.F90
r4911 r5158 251 251 ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we 252 252 ! don't take the logarithm of the first pressure in each column. 253 dojcol = i1,i2253 DO jcol = i1,i2 254 254 this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length) & 255 255 & * thermodynamics%temperature_hl(jcol,2) & … … 258 258 end do 259 259 260 dojlev = 2,nlev-1261 dojcol = i1,i2260 DO jlev = 2,nlev-1 261 DO jcol = i1,i2 262 262 this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) & 263 263 & * thermodynamics%temperature_hl(jcol,jlev+1) & … … 271 271 ! to top-of-atmosphere). In case pressure_hl(:,nlev+1)=0, we 272 272 ! don't take the logarithm of the last pressure in each column. 273 dojlev = 1,nlev-2274 dojcol = i1,i2273 DO jlev = 1,nlev-2 274 DO jcol = i1,i2 275 275 this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length) & 276 276 & * thermodynamics%temperature_hl(jcol,jlev+1) & … … 280 280 end do 281 281 282 dojcol = i1,i2282 DO jcol = i1,i2 283 283 this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length) & 284 284 & * thermodynamics%temperature_hl(jcol,nlev) & … … 340 340 ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we 341 341 ! don't take the logarithm of the first pressure in each column. 342 dojcol = istartcol,iendcol342 DO jcol = istartcol,iendcol 343 343 this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length(jcol)) & 344 344 & * thermodynamics%temperature_hl(jcol,2) & … … 347 347 end do 348 348 349 dojlev = 2,nlev-1350 dojcol = istartcol,iendcol349 DO jlev = 2,nlev-1 350 DO jcol = istartcol,iendcol 351 351 this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol)) & 352 352 & * thermodynamics%temperature_hl(jcol,jlev+1) & … … 360 360 ! to top-of-atmosphere). In case pressure_hl(:,nlev+1)=0, we 361 361 ! don't take the logarithm of the last pressure in each column. 362 dojlev = 1,nlev-2363 dojcol = istartcol,iendcol362 DO jlev = 1,nlev-2 363 DO jcol = istartcol,iendcol 364 364 this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol)) & 365 365 & * thermodynamics%temperature_hl(jcol,jlev+1) & … … 369 369 end do 370 370 371 dojcol = istartcol,iendcol371 DO jcol = istartcol,iendcol 372 372 this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length(jcol)) & 373 373 & * thermodynamics%temperature_hl(jcol,nlev) & … … 423 423 ! top-of-atmosphere to surface). In case pressure_hl(:,1)=0, we 424 424 ! don't take the logarithm of the first pressure in each column. 425 dojcol = istartcol,iendcol425 DO jcol = istartcol,iendcol 426 426 this%overlap_param(jcol,1) = exp(-(R_over_g/decorrelation_length(jcol,1)) & 427 427 & * thermodynamics%temperature_hl(jcol,2) & … … 430 430 end do 431 431 432 dojlev = 2,nlev-1433 dojcol = istartcol,iendcol432 DO jlev = 2,nlev-1 433 DO jcol = istartcol,iendcol 434 434 this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol,jlev)) & 435 435 & * thermodynamics%temperature_hl(jcol,jlev+1) & … … 443 443 ! to top-of-atmosphere). In case pressure_hl(:,nlev+1)=0, we 444 444 ! don't take the logarithm of the last pressure in each column. 445 dojlev = 1,nlev-2446 dojcol = istartcol,iendcol445 DO jlev = 1,nlev-2 446 DO jcol = istartcol,iendcol 447 447 this%overlap_param(jcol,jlev) = exp(-(0.5_jprb*R_over_g/decorrelation_length(jcol,jlev)) & 448 448 & * thermodynamics%temperature_hl(jcol,jlev+1) & … … 452 452 end do 453 453 454 dojcol = istartcol,iendcol454 DO jcol = istartcol,iendcol 455 455 ! AI ATTENTION a verifier decorrelation_length(jcol,nlev-1) ou nlev 456 456 this%overlap_param(jcol,nlev-1) = exp(-(R_over_g/decorrelation_length(jcol,nlev-1)) & … … 663 663 end if 664 664 665 dojcol = i1,i2665 DO jcol = i1,i2 666 666 eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) & 667 667 & * (0.5_jprb / pressure_hl(jcol,isurf)) … … 761 761 end if 762 762 763 dojcol = i1,i2763 DO jcol = i1,i2 764 764 eta = (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)) & 765 765 & * (0.5_jprb / pressure_hl(jcol,isurf)) … … 804 804 ntype = size(this%mixing_ratio,3) 805 805 806 dojlev = 1,nlev807 dojcol = istartcol,iendcol806 DO jlev = 1,nlev 807 DO jcol = istartcol,iendcol 808 808 sum_mixing_ratio(jcol) = 0.0_jprb 809 809 end do 810 dojh = 1, ntype811 dojcol = istartcol,iendcol810 DO jh = 1, ntype 811 DO jcol = istartcol,iendcol 812 812 sum_mixing_ratio(jcol) = sum_mixing_ratio(jcol) + this%mixing_ratio(jcol,jlev,jh) 813 813 end do 814 814 end do 815 dojcol = istartcol,iendcol815 DO jcol = istartcol,iendcol 816 816 if (this%fraction(jcol,jlev) < cloud_fraction_threshold & 817 817 & .or. sum_mixing_ratio(jcol) < cloud_mixing_ratio_threshold) then -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_cover.F90
r4853 r5158 159 159 cum_product = 1.0_jprb - frac(1) 160 160 cum_cloud_cover(1) = frac(1) 161 dojlev = 1,nlev-1161 DO jlev = 1,nlev-1 162 162 ! Compute the combined cloud cover of layers jlev and jlev+1 163 163 pair_cloud_cover(jlev) = max(frac(jlev),frac(jlev+1)) … … 242 242 cum_product = 1.0_jprb - frac(1) 243 243 cum_cloud_cover(1) = frac(1) 244 dojlev = 1,nlev-1244 DO jlev = 1,nlev-1 245 245 ! Convert to "alpha" overlap parameter if necessary 246 246 if (do_overlap_conversion) then … … 258 258 #ifdef DWD_VECTOR_OPTIMIZATIONS 259 259 end do 260 dojlev = 1,nlev-1260 DO jlev = 1,nlev-1 261 261 #endif 262 262 if (frac(jlev) >= MaxCloudFrac) then … … 385 385 jlev = 1 386 386 nobj = 0 387 dowhile (jlev <= nlev)387 DO while (jlev <= nlev) 388 388 if (frac(jlev) > min_frac) then 389 389 ! Starting a new object: store its top … … 392 392 ! Find its maximum cloud fraction 393 393 jlev = jlev + 1 394 do while (jlev <= nlev)394 DO while (jlev <= nlev) 395 395 if (frac(jlev) < frac(jlev-1)) then 396 396 exit … … 400 400 i_max_obj(nobj) = jlev - 1 401 401 ! Find its base 402 dowhile (jlev <= nlev)402 DO while (jlev <= nlev) 403 403 if (frac(jlev) > frac(jlev-1) .or. frac(jlev) <= min_frac) then 404 404 exit … … 435 435 436 436 ! Compute the combined cloud cover of pairs of layers 437 dojlev = 1,nlev-1437 DO jlev = 1,nlev-1 438 438 pair_cloud_cover(jlev) & 439 439 & = overlap_param(jlev)*max(frac(jlev),frac(jlev+1)) & … … 444 444 ! adjacent objects as the product of the layerwise overlap 445 445 ! parameters between their layers of maximum cloud fraction 446 dojobj = 1,nobj-1446 DO jobj = 1,nobj-1 447 447 alpha_obj(jobj) & 448 448 & = product(overlap_param(i_max_obj(jobj):i_max_obj(jobj+1)-1)) … … 456 456 & frac(1:nlev-1), frac(2:nlev)) 457 457 ! Compute the combined cloud cover of pairs of layers 458 dojlev = 1,nlev-1458 DO jlev = 1,nlev-1 459 459 pair_cloud_cover(jlev) & 460 460 & = overlap_alpha(jlev)*max(frac(jlev),frac(jlev+1)) & … … 465 465 ! adjacent objects as the product of the layerwise overlap 466 466 ! parameters between their layers of maximum cloud fraction 467 dojobj = 1,nobj-1467 DO jobj = 1,nobj-1 468 468 alpha_obj(jobj) & 469 469 & = product(overlap_alpha(i_max_obj(jobj):i_max_obj(jobj+1)-1)) … … 476 476 ! of each object: this will later be converted to the cumulative 477 477 ! cloud cover working down from TOA 478 dojobj = 1,nobj478 DO jobj = 1,nobj 479 479 cum_cloud_cover(i_top_obj(jobj)) = frac(i_top_obj(jobj)) 480 dojlev = i_top_obj(jobj), i_base_obj(jobj)-1480 DO jlev = i_top_obj(jobj), i_base_obj(jobj)-1 481 481 if (frac(jlev) >= MaxCloudFrac) then 482 482 ! Cloud cover has reached one … … 496 496 ! Sequentially combine objects until there is only one left 497 497 ! covering the full vertical extent of clouds in the profile 498 dowhile (nobj > 1)498 DO while (nobj > 1) 499 499 ! Find the most correlated adjacent pair of objects 500 500 alpha_max = 0.0_jprb … … 506 506 507 507 jobj = 1 508 dowhile (jobj < nobj)508 DO while (jobj < nobj) 509 509 if (alpha_obj(jobj) > alpha_max) then 510 510 alpha_max = alpha_obj(jobj) … … 533 533 ! Scale the combined cloud cover of the lower object to 534 534 ! account for its overlap with the upper object 535 dojlev = i_top_obj(iobj2),i_base_obj(iobj2)535 DO jlev = i_top_obj(iobj2),i_base_obj(iobj2) 536 536 cum_cloud_cover(jlev) = cum_cloud_cover(i_base_obj(iobj1)) & 537 537 + cum_cloud_cover(jlev) * scaling … … 554 554 ! Ensure that the combined cloud cover of pairs of layers is 555 555 ! consistent with the overhang 556 dojlev = 1,nlev-1556 DO jlev = 1,nlev-1 557 557 pair_cloud_cover(jlev) = max(pair_cloud_cover(jlev), & 558 558 & 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 152 152 153 153 total_cloud_cover = cum_cloud_cover(nlev); 154 dojlev = 1,nlev-1154 DO jlev = 1,nlev-1 155 155 overhang(jlev) = cum_cloud_cover(jlev+1)-cum_cloud_cover(jlev) 156 156 end do … … 166 166 ! Find range of cloudy layers 167 167 jlev = 1 168 do while (frac(jlev) <= 0.0_jprb)168 DO while (frac(jlev) <= 0.0_jprb) 169 169 jlev = jlev + 1 170 170 end do 171 171 ibegin = jlev 172 172 iend = jlev 173 dojlev = jlev+1,nlev173 DO jlev = jlev+1,nlev 174 174 if (frac(jlev) > 0.0_jprb) then 175 175 iend = jlev … … 180 180 overlap_param_inhom = overlap_param 181 181 182 dojlev = ibegin,iend-1182 DO jlev = ibegin,iend-1 183 183 if (overlap_param(jlev) > 0.0_jprb) then 184 184 overlap_param_inhom(jlev) & … … 208 208 209 209 ! Loop over ng columns 210 dojg = 1,ng210 DO jg = 1,ng 211 211 ! Find the cloud top height corresponding to the current 212 212 ! random number, and store in itrigger 213 213 trigger = rand_top(jg) * total_cloud_cover 214 214 jlev = ibegin 215 dowhile (trigger > cum_cloud_cover(jlev) .and. jlev < iend)215 DO while (trigger > cum_cloud_cover(jlev) .and. jlev < iend) 216 216 jlev = jlev + 1 217 217 end do … … 328 328 ! Loop from the layer below the local cloud top down to the 329 329 ! bottom-most cloudy layer 330 dojlev = itrigger+1,iend+1330 DO jlev = itrigger+1,iend+1 331 331 do_fill_od_scaling = .false. 332 332 if (jlev <= iend) then … … 367 367 368 368 ! Loop through the sequence of cloudy layers 369 dojcloud = 2,n_layers_to_scale369 DO jcloud = 2,n_layers_to_scale 370 370 ! Use second random number, and inhomogeneity overlap 371 371 ! parameter, to decide whether the first random number … … 463 463 ! Loop from the layer below the local cloud top down to the 464 464 ! bottom-most cloudy layer 465 dojlev = itrigger+1,iend465 DO jlev = itrigger+1,iend 466 466 iy = iy+1 467 467 if (is_cloudy(jlev-1)) then … … 494 494 495 495 ! Loop through the sequence of cloudy layers 496 dojcloud = 2,n_layers_to_scale496 DO jcloud = 2,n_layers_to_scale 497 497 ! Use second random number, and inhomogeneity overlap 498 498 ! parameter, to decide whether the first random number … … 679 679 680 680 ! Loop down through layers starting at the first cloudy layer 681 dojlev = ibegin,iend681 DO jlev = ibegin,iend 682 682 683 683 if (is_any_cloud(jlev)) then … … 685 685 ! Added for DWD (2020) 686 686 !NEC$ shortloop 687 dojg = 1,ng687 DO jg = 1,ng 688 688 ! The intention is that all these operations are vectorizable, 689 689 ! since all are vector operations on vectors of length ng... -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloud_optics.F90
r4773 r5158 288 288 end if 289 289 290 dojlev = 1,nlev291 dojcol = istartcol,iendcol290 DO jlev = 1,nlev 291 DO jcol = istartcol,iendcol 292 292 ! Only do anything if cloud is present (assume that 293 293 ! cloud%crop_cloud_fraction has already been called) … … 459 459 ! Added for DWD (2020) 460 460 !NEC$ shortloop 461 dojb = 1, config%n_bands_lw461 DO jb = 1, config%n_bands_lw 462 462 od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) + od_lw_ice(jb) 463 463 if (scat_od_lw_liq(jb)+scat_od_lw_ice(jb) > 0.0_jprb) then … … 477 477 ! Added for DWD (2020) 478 478 !NEC$ shortloop 479 dojb = 1, config%n_bands_lw479 DO jb = 1, config%n_bands_lw 480 480 od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) - scat_od_lw_liq(jb) & 481 481 & + od_lw_ice(jb) - scat_od_lw_ice(jb) … … 484 484 ! Added for DWD (2020) 485 485 !NEC$ shortloop 486 dojb = 1, config%n_bands_sw486 DO jb = 1, config%n_bands_sw 487 487 od_sw_cloud(jb,jlev,jcol) = od_sw_liq(jb) + od_sw_ice(jb) 488 488 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 95 95 96 96 ! Loop through columns 97 dojcol = istartcol,iendcol97 DO jcol = istartcol,iendcol 98 98 99 99 ! Compute the reflectance and transmittance of all layers, 100 100 ! neglecting clouds 101 dojlev = 1,nlev101 DO jlev = 1,nlev 102 102 if (config%do_lw_aerosol_scattering) then 103 103 ! Scattering case: first compute clear-sky reflectance, -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_cloudless_sw.F90
r4773 r5158 103 103 104 104 ! Loop through columns 105 dojcol = istartcol,iendcol105 DO jcol = istartcol,iendcol 106 106 ! Only perform calculation if sun above the horizon 107 107 if (single_level%cos_sza(jcol) > 0.0_jprb) then … … 114 114 ! Delta-Eddington scaling has already been performed to the 115 115 ! aerosol part of od, ssa and g 116 dojlev = 1,nlev116 DO jlev = 1,nlev 117 117 call calc_two_stream_gammas_sw(ng, cos_sza, & 118 118 & ssa(:,jlev,jcol), g(:,jlev,jcol), & … … 128 128 else 129 129 ! Apply delta-Eddington scaling to the aerosol-gas mixture 130 dojlev = 1,nlev130 DO jlev = 1,nlev 131 131 od_total = od(:,jlev,jcol) 132 132 ssa_total = ssa(:,jlev,jcol) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_config.F90
r4853 r5158 1657 1657 ! contains the weights of interest. We now find the non-zero weights 1658 1658 nweights = 0 1659 dojband = 1,size(mapping,2)1659 DO jband = 1,size(mapping,2) 1660 1660 if (mapping(2,jband) > 0.0_jprb) then 1661 1661 nweights = nweights+1 … … 1674 1674 & wavenumber2, ' cm-1):' 1675 1675 if (this%do_cloud_aerosol_per_sw_g_point) then 1676 dojband = 1, nweights1676 DO jband = 1, nweights 1677 1677 write(nulout, '(a,i0,a,f8.4)') ' Shortwave g point ', iband(jband), ': ', weight(jband) 1678 1678 end do 1679 1679 else 1680 dojband = 1, nweights1680 DO jband = 1, nweights 1681 1681 wavenumber1_band = this%gas_optics_sw%spectral_def%wavenumber1_band(iband(jband)) 1682 1682 wavenumber2_band = this%gas_optics_sw%spectral_def%wavenumber2_band(iband(jband)) … … 1728 1728 allocate(diag_ind(ninterval+2)) 1729 1729 diag_ind = 0 1730 dojint = 1,ninterval+21730 DO jint = 1,ninterval+2 1731 1731 diag_ind(jint) = jint 1732 1732 end do … … 1894 1894 ! Count the number of albedo/emissivity intervals 1895 1895 ninterval = 0 1896 dojint = 1,NMaxAlbedoIntervals1896 DO jint = 1,NMaxAlbedoIntervals 1897 1897 if (this%i_sw_albedo_index(jint) > 0) then 1898 1898 ninterval = jint … … 1948 1948 write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_albedo_from_band_sw), & 1949 1949 & ' shortwave intervals to albedo intervals:' 1950 dojband = 1,size(this%i_albedo_from_band_sw)1950 DO jband = 1,size(this%i_albedo_from_band_sw) 1951 1951 write(nulout,'(a,i0)',advance='no') ' ', this%i_albedo_from_band_sw(jband) 1952 1952 end do … … 1972 1972 ! Count the number of albedo/emissivity intervals 1973 1973 ninterval = 0 1974 dojint = 1,NMaxAlbedoIntervals1974 DO jint = 1,NMaxAlbedoIntervals 1975 1975 if (this%i_lw_emiss_index(jint) > 0) then 1976 1976 ninterval = jint … … 2026 2026 write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_emiss_from_band_lw), & 2027 2027 & ' longwave intervals to emissivity intervals:' 2028 dojband = 1,size(this%i_emiss_from_band_lw)2028 DO jband = 1,size(this%i_emiss_from_band_lw) 2029 2029 write(nulout,'(a,i0)',advance='no') ' ', this%i_emiss_from_band_lw(jband) 2030 2030 end do … … 2056 2056 is_not_found = .true. 2057 2057 2058 dojc = 0,size(enum_str)-12058 DO jc = 0,size(enum_str)-1 2059 2059 if (trim(str) == trim(enum_str(jc))) then 2060 2060 icode = jc … … 2066 2066 write(nulerr,'(a,a,a,a,a)',advance='no') '*** Error: ', trim(var_name), & 2067 2067 & ' must be one of: "', enum_str(0), '"' 2068 dojc = 1,size(enum_str)-12068 DO jc = 1,size(enum_str)-1 2069 2069 write(nulerr,'(a,a,a)',advance='no') ', "', trim(enum_str(jc)), '"' 2070 2070 end do -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_delta_eddington.h
r4773 r5158 86 86 integer :: j 87 87 88 doj = 1,ng88 DO j = 1,ng 89 89 g = scat_od_g(j) / max(scat_od(j), 1.0e-24) 90 90 f = g*g -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd.F90
r4853 r5158 209 209 istart = 1 210 210 this%i_gas_mapping = 0 211 dojgas = 1, this%ngas211 DO jgas = 1, this%ngas 212 212 if (jgas < this%ngas) then 213 213 inext = istart + scan(constituent_id(istart:nchar), ' ') … … 217 217 ! Find gas code 218 218 i_gas_code = 0 219 dojjgas = 1, NMaxGases219 DO jjgas = 1, NMaxGases 220 220 if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then 221 221 i_gas_code = jjgas … … 264 264 & this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries' 265 265 write(nulout, '(a)') ' Gases:' 266 dojgas = 1,this%ngas266 DO jgas = 1,this%ngas 267 267 if (this%single_gas(jgas)%i_gas_code > 0) then 268 268 write(nulout, '(a,a,a)', advance='no') ' ', & … … 373 373 374 374 ! Interpolate input SSI to regular wavenumber grid 375 dojwav_grid = 1,nwav_grid376 dojwav = 1,nwav-1375 DO jwav_grid = 1,nwav_grid 376 DO jwav = 1,nwav-1 377 377 if (wavenumber(jwav) < wavenumber_grid(jwav_grid) & 378 378 & .and. wavenumber(jwav+1) >= wavenumber_grid(jwav_grid)) then … … 426 426 & ReferenceTSI, ' W m-2, solar cycle amplitude (at solar maximum), update to original solar irradiance' 427 427 iband = 0 428 dojg = 1,this%ng428 DO jg = 1,this%ng 429 429 if (this%spectral_def%i_band_number(jg) > iband) then 430 430 iband = this%spectral_def%i_band_number(jg) … … 516 516 global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass) 517 517 518 dojcol = istartcol,iendcol518 DO jcol = istartcol,iendcol 519 519 520 520 log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1))) 521 521 522 dojlev = 1,nlev522 DO jlev = 1,nlev 523 523 ! Find interpolation points in pressure 524 524 pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) & … … 546 546 optical_depth_fl(:,:,jcol) = 0.0_jprb 547 547 548 dojgas = 1,this%ngas548 DO jgas = 1,this%ngas 549 549 550 550 associate (single_gas => this%single_gas(jgas)) … … 561 561 end if 562 562 563 dojlev = 1,nlev563 DO jlev = 1,nlev 564 564 optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & 565 565 & + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & … … 581 581 end if 582 582 583 dojlev = 1,nlev583 DO jlev = 1,nlev 584 584 optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & 585 585 & + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & … … 592 592 ! Composite gases 593 593 molar_abs => this%single_gas(jgas)%molar_abs 594 dojlev = 1,nlev594 DO jlev = 1,nlev 595 595 optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & 596 596 & + (simple_multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & … … 611 611 molar_abs_conc => this%single_gas(jgas)%molar_abs_conc 612 612 mole_frac1 = exp(single_gas%log_mole_frac1) 613 dojlev = 1,nlev613 DO jlev = 1,nlev 614 614 ! Take care of mole_fraction == 0 615 615 log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode)*scaling, mole_frac1)) … … 628 628 ! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) & 629 629 ! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1))) 630 dojlev = 1,nlev630 DO jlev = 1,nlev 631 631 optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & 632 632 & + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode) * scaling) * ( & … … 651 651 ! Rayleigh scattering 652 652 if (this%is_sw .and. present(rayleigh_od_fl)) then 653 dojlev = 1,nlev653 DO jlev = 1,nlev 654 654 rayleigh_od_fl(:,jlev,jcol) = global_multiplier & 655 655 & * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat … … 729 729 od_fl(:,:,:) = 0.0_jprb 730 730 731 dojlev = 1,nlev732 dojcol = istartcol,iendcol731 DO jlev = 1,nlev 732 DO jcol = istartcol,iendcol 733 733 734 734 log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,jlev)+pressure_hl(jcol,jlev+1))) … … 758 758 end do 759 759 760 dojgas = 1,this%ngas760 DO jgas = 1,this%ngas 761 761 762 762 associate (single_gas => this%single_gas(jgas)) … … 768 768 molar_abs => this%single_gas(jgas)%molar_abs 769 769 770 dojlev = 1,nlev771 dojg = 1, this%ng772 dojcol = istartcol,iendcol770 DO jlev = 1,nlev 771 DO jg = 1, this%ng 772 DO jcol = istartcol,iendcol 773 773 multiplier = simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode) 774 774 … … 785 785 molar_abs => this%single_gas(jgas)%molar_abs 786 786 787 dojlev = 1,nlev788 dojg = 1, this%ng789 dojcol = istartcol,iendcol787 DO jlev = 1,nlev 788 DO jg = 1, this%ng 789 DO jcol = istartcol,iendcol 790 790 multiplier = simple_multiplier(jcol,jlev) * (mole_fraction_fl(jcol,jlev,igascode) & 791 791 & - single_gas%reference_mole_frac) … … 804 804 molar_abs => this%single_gas(jgas)%molar_abs 805 805 806 dojlev = 1,nlev807 dojg = 1, this%ng808 dojcol = istartcol,iendcol806 DO jlev = 1,nlev 807 DO jg = 1, this%ng 808 DO jcol = istartcol,iendcol 809 809 od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) & 810 810 & + (simple_multiplier(jcol,jlev)*tw1(jcol,jlev)) * & … … 823 823 mole_frac1 = exp(single_gas%log_mole_frac1) 824 824 825 dojlev = 1,nlev826 dojcol = istartcol,iendcol825 DO jlev = 1,nlev 826 DO jcol = istartcol,iendcol 827 827 ! Take care of mole_fraction == 0 828 828 log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1)) … … 835 835 end do 836 836 837 dojlev = 1,nlev838 dojg = 1, this%ng837 DO jlev = 1,nlev 838 DO jg = 1, this%ng 839 839 !NEC$ select_vector 840 dojcol = istartcol,iendcol840 DO jcol = istartcol,iendcol 841 841 842 842 od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) & … … 866 866 867 867 ! Ensure the optical depth is not negative 868 dojcol = istartcol,iendcol869 dojlev = 1,nlev870 dojg = 1, this%ng868 DO jcol = istartcol,iendcol 869 DO jlev = 1,nlev 870 DO jg = 1, this%ng 871 871 optical_depth_fl(jg,jlev,jcol) = max(0.0_jprb, od_fl(jcol,jg,jlev)) 872 872 end do … … 876 876 ! Rayleigh scattering 877 877 if (this%is_sw .and. present(rayleigh_od_fl)) then 878 dojcol = istartcol,iendcol879 dojlev = 1,nlev880 dojg = 1, this%ng878 DO jcol = istartcol,iendcol 879 DO jlev = 1,nlev 880 DO jg = 1, this%ng 881 881 rayleigh_od_fl(jg,jlev,jcol) = global_multiplier & 882 882 & * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat(jg) … … 908 908 integer :: it1, jt 909 909 910 dojt = 1,nt910 DO jt = 1,nt 911 911 tindex1 = (temperature(jt) - this%temperature1_planck) & 912 912 & * (1.0_jprb / this%d_temperature_planck) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ecckd_interface.F90
r4853 r5158 272 272 ! rayleigh optical depth: convert to total optical depth and 273 273 ! single-scattering albedo 274 dojcol = istartcol,iendcol275 dojlev = 1, nlev276 dojg = 1, config%n_g_sw274 DO jcol = istartcol,iendcol 275 DO jlev = 1, nlev 276 DO jg = 1, config%n_g_sw 277 277 od_sw(jg,jlev,jcol) = od_sw(jg,jlev,jcol) + ssa_sw(jg,jlev,jcol) 278 278 ssa_sw(jg,jlev,jcol) = ssa_sw(jg,jlev,jcol) / od_sw(jg,jlev,jcol) … … 308 308 309 309 ! Calculate the Planck function for each g point 310 dojcol = istartcol,iendcol310 DO jcol = istartcol,iendcol 311 311 call config%gas_optics_lw%calc_planck_function(nlev+1, & 312 312 & thermodynamics%temperature_hl(jcol,:), planck_hl(:,:,jcol)) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_flux.F90
r4853 r5158 423 423 & config%i_band_from_reordered_g_sw, & 424 424 & this%sw_dn_surf_band, istartcol, iendcol) 425 dojcol = istartcol,iendcol425 DO jcol = istartcol,iendcol 426 426 this%sw_dn_surf_band(:,jcol) & 427 427 & = this%sw_dn_surf_band(:,jcol) & … … 429 429 end do 430 430 else 431 dojcol = istartcol,iendcol431 DO jcol = istartcol,iendcol 432 432 call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), & 433 433 & config%i_band_from_reordered_g_sw, & … … 450 450 & config%i_band_from_reordered_g_sw, & 451 451 & this%sw_dn_surf_clear_band, istartcol, iendcol) 452 dojcol = istartcol,iendcol452 DO jcol = istartcol,iendcol 453 453 this%sw_dn_surf_clear_band(:,jcol) & 454 454 & = this%sw_dn_surf_clear_band(:,jcol) & … … 456 456 end do 457 457 else 458 dojcol = istartcol,iendcol458 DO jcol = istartcol,iendcol 459 459 call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), & 460 460 & config%i_band_from_reordered_g_sw, & … … 486 486 & this%sw_dn_diffuse_surf_canopy, istartcol, iendcol) 487 487 else 488 dojcol = istartcol,iendcol488 DO jcol = istartcol,iendcol 489 489 call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), & 490 490 & config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), & … … 502 502 this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = 0.0_jprb 503 503 this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = 0.0_jprb 504 dojband = 1,config%n_bands_sw505 dojalbedoband = 1,nalbedoband504 DO jband = 1,config%n_bands_sw 505 DO jalbedoband = 1,nalbedoband 506 506 if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then 507 507 ! Initially, "diffuse" is actually "total" … … 534 534 & this%lw_dn_surf_canopy, istartcol, iendcol) 535 535 else 536 dojcol = istartcol,iendcol536 DO jcol = istartcol,iendcol 537 537 call indexed_sum(this%lw_dn_surf_g(:,jcol), & 538 538 & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), & … … 548 548 & lw_dn_surf_band, istartcol, iendcol) 549 549 else 550 dojcol = istartcol,iendcol550 DO jcol = istartcol,iendcol 551 551 call indexed_sum(this%lw_dn_surf_g(:,jcol), & 552 552 & config%i_band_from_reordered_g_lw, & … … 556 556 nalbedoband = size(config%lw_emiss_weights,1) 557 557 this%lw_dn_surf_canopy(:,istartcol:iendcol) = 0.0_jprb 558 dojband = 1,config%n_bands_lw559 dojalbedoband = 1,nalbedoband558 DO jband = 1,config%n_bands_lw 559 DO jalbedoband = 1,nalbedoband 560 560 if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then 561 561 this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) & … … 602 602 & this%sw_up_toa_band, istartcol, iendcol) 603 603 else 604 dojcol = istartcol,iendcol604 DO jcol = istartcol,iendcol 605 605 call indexed_sum(this%sw_dn_toa_g(:,jcol), & 606 606 & config%i_band_from_reordered_g_sw, & … … 618 618 & this%sw_up_toa_clear_band, istartcol, iendcol) 619 619 else 620 dojcol = istartcol,iendcol620 DO jcol = istartcol,iendcol 621 621 call indexed_sum(this%sw_up_toa_clear_g(:,jcol), & 622 622 & config%i_band_from_reordered_g_sw, & … … 634 634 & this%lw_up_toa_band, istartcol, iendcol) 635 635 else 636 dojcol = istartcol,iendcol636 DO jcol = istartcol,iendcol 637 637 call indexed_sum(this%lw_up_toa_g(:,jcol), & 638 638 & config%i_band_from_reordered_g_lw, & … … 647 647 & this%lw_up_toa_clear_band, istartcol, iendcol) 648 648 else 649 dojcol = istartcol,iendcol649 DO jcol = istartcol,iendcol 650 650 call indexed_sum(this%lw_up_toa_clear_g(:,jcol), & 651 651 & config%i_band_from_reordered_g_lw, & … … 753 753 iend = ubound(source,1) 754 754 755 dojg = istart, iend755 DO jg = istart, iend 756 756 ig = ind(jg) 757 757 dest(ig) = dest(ig) + source(jg) … … 777 777 iend = ubound(source,1) 778 778 779 dojg = istart, iend779 DO jg = istart, iend 780 780 ig = ind(jg) 781 781 dest(ig) = dest(ig) + source(jg) … … 797 797 dest = 0.0 798 798 799 dojg = lbound(source,1), ubound(source,1)799 DO jg = lbound(source,1), ubound(source,1) 800 800 ig = ind(jg) 801 dojc = ist, iend801 DO jc = ist, iend 802 802 dest(ig,jc) = dest(ig,jc) + source(jg,jc) 803 803 end do … … 820 820 nlev = size(source,2) 821 821 822 dojlev = 1,nlev823 dojg = istart, iend822 DO jlev = 1,nlev 823 DO jg = istart, iend 824 824 ig = ind(jg) 825 825 dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev) … … 846 846 nlev = size(source,2) 847 847 848 dojlev = 1,nlev849 dojg = istart, iend848 DO jlev = 1,nlev 849 DO jg = istart, iend 850 850 ig = ind(jg) 851 851 dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_gas.F90
r4853 r5158 209 209 this%is_well_mixed(igas) = .false. 210 210 211 dojk = 1,this%nlev212 dojc = i1,i2211 DO jk = 1,this%nlev 212 DO jc = i1,i2 213 213 this%mixing_ratio(jc,jk,igas) = mixing_ratio(jc-i1+1,jk) 214 214 end do … … 294 294 this%is_well_mixed(igas) = .true. 295 295 296 dojk = 1,this%nlev297 dojc = i1,i2296 DO jk = 1,this%nlev 297 DO jc = i1,i2 298 298 this%mixing_ratio(jc,jk,igas) = mixing_ratio 299 299 end do … … 397 397 end if 398 398 else 399 dojg = 1,this%ntype399 DO jg = 1,this%ntype 400 400 call this%set_units(iunits, igas=this%icode(jg), scale_factor=new_sf) 401 401 end do … … 416 416 417 417 scaling = this%scale_factor 418 dojg = 1,NMaxGases418 DO jg = 1,NMaxGases 419 419 if (iunits == IMassMixingRatio .and. this%iunits(jg) == IVolumeMixingRatio) then 420 420 scaling(jg) = scaling(jg) * GasMolarMass(jg) / AirMolarMass … … 481 481 end if 482 482 else 483 dojg = 1,this%ntype483 DO jg = 1,this%ntype 484 484 call this%assert_units(iunits, igas=this%icode(jg), scale_factor=sf, istatus=istatus) 485 485 end do -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics.F90
r4853 r5158 53 53 ! Count number of cloud types 54 54 config%n_cloud_types = 0 55 dojtype = 1,NMaxCloudTypes55 DO jtype = 1,NMaxCloudTypes 56 56 if (len_trim(config%cloud_type_name(jtype)) > 0) then 57 57 config%n_cloud_types = jtype … … 83 83 84 84 ! Load cloud optics data 85 dojtype = 1,config%n_cloud_types85 DO jtype = 1,config%n_cloud_types 86 86 if (config%cloud_type_name(jtype)(1:1) == '/') then 87 87 file_name = trim(config%cloud_type_name(jtype)) … … 191 191 192 192 ! Loop over cloud types 193 dojtype = 1,config%n_cloud_types193 DO jtype = 1,config%n_cloud_types 194 194 ! Compute in-cloud water path 195 195 if (config%is_homogeneous) then … … 238 238 ! Scale the combined longwave optical properties 239 239 if (config%do_lw_cloud_scattering) then 240 dojcol = istartcol, iendcol241 dojlev = 1,nlev240 DO jcol = istartcol, iendcol 241 DO jlev = 1,nlev 242 242 if (cloud%fraction(jcol,jlev) > 0.0_jprb) then 243 243 ! Note that original cloud optics does not do … … 259 259 if (config%do_sw) then 260 260 if (.not. config%do_sw_delta_scaling_with_gases) then 261 dojcol = istartcol, iendcol262 dojlev = 1,nlev261 DO jcol = istartcol, iendcol 262 DO jlev = 1,nlev 263 263 if (cloud%fraction(jcol,jlev) > 0.0_jprb) then 264 264 call delta_eddington_extensive(od_sw_cloud(:,jlev,jcol), & … … 269 269 end if 270 270 271 dojcol = istartcol, iendcol272 dojlev = 1,nlev271 DO jcol = istartcol, iendcol 272 DO jlev = 1,nlev 273 273 if (cloud%fraction(jcol,jlev) > 0.0_jprb) then 274 274 ! Scale to get asymmetry factor and single scattering albedo 275 dojg = 1, config%n_bands_sw275 DO jg = 1, config%n_bands_sw 276 276 g_sw_cloud(jg,jlev,jcol) = g_sw_cloud(jg,jlev,jcol) & 277 277 & / max(ssa_sw_cloud(jg,jlev,jcol), 1.0e-15_jprb) … … 308 308 if (lhook) call dr_hook('radiation_general_cloud_optics:save',0,hook_handle) 309 309 310 dojtype = 1,config%n_cloud_types310 DO jtype = 1,config%n_cloud_types 311 311 if (config%do_sw) then 312 312 associate(co_sw => config%cloud_optics_sw(jtype)) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics_data.F90
r4936 r5158 285 285 286 286 if (present(scat_od)) then 287 dojcol = 1,ncol288 dojlev = 1,nlev287 DO jcol = 1,ncol 288 DO jlev = 1,nlev 289 289 if (cloud_fraction(jcol, jlev) > 0.0_jprb) then 290 290 re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) & … … 293 293 weight2 = re_index - ire 294 294 weight1 = 1.0_jprb - weight2 295 dojg = 1, ng295 DO jg = 1, ng 296 296 od_local = water_path(jcol, jlev) * (weight1*this%mass_ext(jg,ire) & 297 297 & +weight2*this%mass_ext(jg,ire+1)) … … 310 310 else 311 311 ! No scattering: return the absorption optical depth 312 dojcol = 1,ncol313 dojlev = 1,nlev312 DO jcol = 1,ncol 313 DO jlev = 1,nlev 314 314 if (water_path(jcol, jlev) > 0.0_jprb) then 315 315 re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) & … … 404 404 405 405 ! Define effective radius 406 doire = 1,this%n_effective_radius406 DO ire = 1,this%n_effective_radius 407 407 effective_radius(ire) = this%effective_radius_0 + this%d_effective_radius*(ire-1) 408 408 end do -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_lw.F90
r4773 r5158 119 119 120 120 ! Loop through columns 121 dojcol = istartcol,iendcol121 DO jcol = istartcol,iendcol 122 122 123 123 ! Is there any cloud in the profile? 124 124 is_cloudy_profile = .false. 125 dojlev = 1,nlev125 DO jlev = 1,nlev 126 126 if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then 127 127 is_cloudy_profile = .true. … … 135 135 ! the clear-sky layers since these will be needed when we come 136 136 ! to do the total-sky fluxes. 137 dojlev = 1,nlev137 DO jlev = 1,nlev 138 138 if (config%do_clear .or. cloud%fraction(jcol,jlev) & 139 139 & < config%cloud_fraction_threshold) then … … 201 201 ! now. 202 202 if (is_cloudy_profile .or. .not. config%do_clear) then 203 dojlev = 1,nlev203 DO jlev = 1,nlev 204 204 ! Compute combined gas+aerosol+cloud optical properties; 205 205 ! note that for clear layers, the reflectance and -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_homogeneous_sw.F90
r4773 r5158 123 123 124 124 ! Loop through columns 125 dojcol = istartcol,iendcol125 DO jcol = istartcol,iendcol 126 126 ! Only perform calculation if sun above the horizon 127 127 if (single_level%cos_sza(jcol) > 0.0_jprb) then … … 131 131 ! Is there any cloud in the profile? 132 132 is_cloudy_profile = .false. 133 dojlev = 1,nlev133 DO jlev = 1,nlev 134 134 if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then 135 135 is_cloudy_profile = .true. … … 146 146 ! Delta-Eddington scaling has already been performed to the 147 147 ! aerosol part of od, ssa and g 148 dojlev = 1,nlev148 DO jlev = 1,nlev 149 149 if (config%do_clear .or. cloud%fraction(jcol,jlev) & 150 150 & < config%cloud_fraction_threshold) then … … 164 164 else 165 165 ! Apply delta-Eddington scaling to the aerosol-gas mixture 166 dojlev = 1,nlev166 DO jlev = 1,nlev 167 167 if (config%do_clear .or. cloud%fraction(jcol,jlev) & 168 168 & < config%cloud_fraction_threshold) then … … 229 229 ! fluxes now. 230 230 if (is_cloudy_profile .or. .not. config%do_clear) then 231 dojlev = 1,nlev231 DO jlev = 1,nlev 232 232 ! Compute combined gas+aerosol+cloud optical properties; 233 233 ! note that for clear layers, the reflectance and -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ice_optics_fu.F90
r4773 r5158 73 73 ! Added for DWD (2020) 74 74 !NEC$ shortloop 75 dojb = 1, nb75 DO jb = 1, nb 76 76 od(jb) = iwp_gm_2 * (coeff(jb,1) + coeff(jb,2) * inv_de_um) 77 77 scat_od(jb) = od(jb) * (1.0_jprb - (coeff(jb,3) + de_um*(coeff(jb,4) & … … 124 124 ! Added for DWD (2020) 125 125 !NEC$ shortloop 126 dojb = 1, nb126 DO jb = 1, nb 127 127 od(jb) = iwp_gm_2 * (coeff(jb,1) + inv_de_um*(coeff(jb,2) & 128 128 & + inv_de_um*coeff(jb,3))) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_ifs_rrtm.F90
r4853 r5158 386 386 ZONEMINUS_ARRAY = ZONEMINUS 387 387 388 dojlev=1,nlev389 dojcol= istartcol,iendcol388 DO jlev=1,nlev 389 DO jcol= istartcol,iendcol 390 390 pressure_fl(jcol,jlev) & 391 391 & = 0.5_jprb * (thermodynamics%pressure_hl(jcol,jlev+istartlev-1) & … … 486 486 ! arrays have dimensions in a different order to the inputs, 487 487 ! so there is some inefficiency here. 488 dojgreorder = 1,config%n_g_lw488 DO jgreorder = 1,config%n_g_lw 489 489 iband = config%i_band_from_reordered_g_lw(jgreorder) 490 490 ig = config%i_g_from_reordered_g_lw(jgreorder) 491 491 492 492 ! Top-of-atmosphere half level 493 dojlev = 1,nlev494 dojcol = istartcol,iendcol493 DO jlev = 1,nlev 494 DO jcol = istartcol,iendcol 495 495 ! Some g points can return negative optical depths; 496 496 ! specifically original g points 54-56 which causes … … 504 504 else 505 505 ! G points have not been reordered 506 dojcol = istartcol,iendcol507 dojlev = 1,nlev506 DO jcol = istartcol,iendcol 507 DO jlev = 1,nlev 508 508 ! Check for negative optical depth 509 509 od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol)) … … 544 544 ! Scale the incoming solar per band, if requested 545 545 if (config%use_spectral_solar_scaling) then 546 dojg = 1,JPGPT_SW547 do jcol = istartcol,iendcol546 DO jg = 1,JPGPT_SW 547 DO jcol = istartcol,iendcol 548 548 ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * & 549 549 & single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg)) … … 556 556 ! ZINCSOL will be zero. 557 557 if (present(incoming_sw)) then 558 dojcol = istartcol,iendcol558 DO jcol = istartcol,iendcol 559 559 if (single_level%cos_sza(jcol) > 0.0_jprb) then 560 560 ! Added for DWD (2020) … … 570 570 ! if (.true.) then 571 571 ! Account for reordered g points 572 dojgreorder = 1,config%n_g_sw572 DO jgreorder = 1,config%n_g_sw 573 573 ig = config%i_g_from_reordered_g_sw(jgreorder) 574 dojlev = 1,nlev575 dojcol = istartcol,iendcol574 DO jlev = 1,nlev 575 DO jcol = istartcol,iendcol 576 576 ! Check for negative optical depth 577 577 od_sw (jgreorder,nlev+1-jlev,jcol) & … … 587 587 else 588 588 ! G points have not been reordered 589 dojcol = istartcol,iendcol590 dojlev = 1,nlev591 dojg = 1,config%n_g_sw589 DO jcol = istartcol,iendcol 590 DO jlev = 1,nlev 591 DO jg = 1,config%n_g_sw 592 592 ! Check for negative optical depth 593 593 od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg)) … … 596 596 end do 597 597 if (present(incoming_sw)) then 598 dojg = 1,config%n_g_sw598 DO jg = 1,config%n_g_sw 599 599 incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg) 600 600 end do … … 667 667 ! lowest interpolation bound, and the fraction into interpolation 668 668 ! interval 669 dojlev = 1,nlev+1670 dojcol = istartcol,iendcol669 DO jlev = 1,nlev+1 670 DO jcol = istartcol,iendcol 671 671 temperature = thermodynamics%temperature_hl(jcol,jlev+ilevoffset) 672 672 if (temperature < 339.0_jprb .and. temperature >= 160.0_jprb) then … … 687 687 688 688 ! Calculate Planck functions per band 689 dojband = 1,config%n_bands_lw689 DO jband = 1,config%n_bands_lw 690 690 factor = zfluxfac * delwave(jband) 691 dojcol = istartcol,iendcol691 DO jcol = istartcol,iendcol 692 692 planck_store(jcol,jband) = factor & 693 693 & * (totplnk(ind(jcol),jband) & … … 706 706 ! Top-of-atmosphere half level - note that PFRAC is on model 707 707 ! levels not half levels 708 dojgreorder = 1,config%n_g_lw708 DO jgreorder = 1,config%n_g_lw 709 709 iband = config%i_band_from_reordered_g_lw(jgreorder) 710 710 ig = config%i_g_from_reordered_g_lw(jgreorder) … … 713 713 end do 714 714 else 715 dojgreorder = 1,config%n_g_lw715 DO jgreorder = 1,config%n_g_lw 716 716 iband = config%i_band_from_reordered_g_lw(jgreorder) 717 717 ig = config%i_g_from_reordered_g_lw(jgreorder) … … 726 726 ! Top-of-atmosphere half level - note that PFRAC is on model 727 727 ! levels not half levels 728 dojg = 1,config%n_g_lw728 DO jg = 1,config%n_g_lw 729 729 iband = config%i_band_from_g_lw(jg) 730 730 planck_hl(jg,1,:) = planck_store(:,iband) * PFRAC(:,jg,nlev) 731 731 end do 732 732 else 733 dojg = 1,config%n_g_lw733 DO jg = 1,config%n_g_lw 734 734 iband = config%i_band_from_g_lw(jg) 735 735 planck_tmp(:,jg) = planck_store(:,iband) * PFRAC(:,jg,nlev+2-jlev) 736 736 end do 737 dojcol = istartcol,iendcol737 DO jcol = istartcol,iendcol 738 738 planck_hl(:,jlev,jcol) = planck_tmp(jcol,:) 739 739 end do … … 794 794 795 795 ! Work out surface interpolations 796 dojcol = istartcol,iendcol796 DO jcol = istartcol,iendcol 797 797 Tsurf = temperature(jcol) 798 798 if (Tsurf < 339.0_jprb .and. Tsurf >= 160.0_jprb) then … … 813 813 814 814 ! Calculate Planck functions per band 815 dojband = 1,config%n_bands_lw815 DO jband = 1,config%n_bands_lw 816 816 factor = zfluxfac * delwave(jband) 817 dojcol = istartcol,iendcol817 DO jcol = istartcol,iendcol 818 818 planck_store(jcol,jband) = factor & 819 819 & * (totplnk(ind(jcol),jband) & … … 829 829 ! in pressure since the the functions above treat pressure 830 830 ! decreasing with increasing index. 831 dojgreorder = 1,config%n_g_lw831 DO jgreorder = 1,config%n_g_lw 832 832 iband = config%i_band_from_reordered_g_lw(jgreorder) 833 833 ig = config%i_g_from_reordered_g_lw(jgreorder) … … 836 836 else 837 837 ! G points have not been reordered 838 dojg = 1,config%n_g_lw838 DO jg = 1,config%n_g_lw 839 839 iband = config%i_band_from_g_lw(jg) 840 840 planck_surf(jg,:) = planck_store(:,iband) * PFRAC(:,jg) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_liquid_optics_socrates.F90
r4773 r5158 65 65 ! Added for DWD (2020) 66 66 !NEC$ shortloop 67 dojb = 1, nb67 DO jb = 1, nb 68 68 od(jb) = lwp * (coeff(jb,1) + re*(coeff(jb,2) + re*coeff(jb,3))) & 69 69 & / (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 74 74 ! Move up through the atmosphere computing the derivatives at each 75 75 ! half-level 76 dojlev = nlev,1,-176 DO jlev = nlev,1,-1 77 77 lw_derivatives_g = lw_derivatives_g * transmittance(:,jlev) 78 78 lw_derivatives(icol,jlev) = sum(lw_derivatives_g) … … 122 122 ! Move up through the atmosphere computing the derivatives at each 123 123 ! half-level 124 dojlev = nlev,1,-1124 DO jlev = nlev,1,-1 125 125 lw_derivatives_g = lw_derivatives_g * transmittance(:,jlev) 126 126 lw_derivatives(icol,jlev) = (1.0_jprb - weight) * lw_derivatives(icol,jlev) & … … 177 177 ! Move up through the atmosphere computing the derivatives at each 178 178 ! half-level 179 dojlev = nlev,1,-1179 DO jlev = nlev,1,-1 180 180 ! Compute effect of overlap at half-level jlev+1, yielding 181 181 ! derivatives just above that half-level … … 243 243 ! multiplication inline 244 244 245 dojlev = nlev,1,-1245 DO jlev = nlev,1,-1 246 246 ! Compute effect of overlap at half-level jlev+1, yielding 247 247 ! derivatives just above that half-level … … 249 249 250 250 associate(A=>u_matrix(:,:,jlev+1), b=>lw_deriv_below) 251 do jg = 1,ng251 DO jg = 1,ng 252 252 ! Both inner and outer loop of the matrix loops j1 and j2 unrolled 253 253 ! inner loop: j2=1 j2=2 j2=3 … … 273 273 ! Move up through the atmosphere computing the derivatives at each 274 274 ! half-level 275 dojlev = nlev,1,-1275 DO jlev = nlev,1,-1 276 276 ! Compute effect of overlap at half-level jlev+1, yielding 277 277 ! derivatives just above that half-level -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_matrix.F90
r4773 r5158 91 91 mat_x_vec(1:iend,1) = A(1:iend,1,1)*b(1:iend,1) 92 92 else 93 doj1 = 1,m94 doj2 = 1,m93 DO j1 = 1,m 94 DO j2 = 1,m 95 95 mat_x_vec(1:iend,j1) = mat_x_vec(1:iend,j1) & 96 96 & + A(1:iend,j1,j2)*b(1:iend,j2) … … 125 125 singlemat_x_vec = 0.0_jprb 126 126 127 doj1 = 1,m128 doj2 = 1,m127 DO j1 = 1,m 128 DO j2 = 1,m 129 129 singlemat_x_vec(1:iend,j1) = singlemat_x_vec(1:iend,j1) & 130 130 & + A(j1,j2)*b(1:iend,j2) … … 176 176 m2block = 2*mblock 177 177 ! Do the top-left (C, D, F, G) 178 doj2 = 1,m2block179 doj1 = 1,m2block180 doj3 = 1,m2block178 DO j2 = 1,m2block 179 DO j1 = 1,m2block 180 DO j3 = 1,m2block 181 181 mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) & 182 182 & + A(1:iend,j1,j3)*B(1:iend,j3,j2) … … 184 184 end do 185 185 end do 186 doj2 = m2block+1,m186 DO j2 = m2block+1,m 187 187 ! Do the top-right (E & H) 188 doj1 = 1,m2block189 doj3 = 1,m188 DO j1 = 1,m2block 189 DO j3 = 1,m 190 190 mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) & 191 191 & + A(1:iend,j1,j3)*B(1:iend,j3,j2) … … 193 193 end do 194 194 ! Do the bottom-right (I) 195 doj1 = m2block+1,m196 doj3 = m2block+1,m195 DO j1 = m2block+1,m 196 DO j3 = m2block+1,m 197 197 mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) & 198 198 & + A(1:iend,j1,j3)*B(1:iend,j3,j2) … … 202 202 else 203 203 ! Ordinary dense matrix 204 doj2 = 1,m205 doj1 = 1,m206 doj3 = 1,m204 DO j2 = 1,m 205 DO j1 = 1,m 206 DO j3 = 1,m 207 207 mat_x_mat(1:iend,j1,j2) = mat_x_mat(1:iend,j1,j2) & 208 208 & + A(1:iend,j1,j3)*B(1:iend,j3,j2) … … 238 238 singlemat_x_mat = 0.0_jprb 239 239 240 doj2 = 1,m241 doj1 = 1,m242 doj3 = 1,m240 DO j2 = 1,m 241 DO j1 = 1,m 242 DO j3 = 1,m 243 243 singlemat_x_mat(1:iend,j1,j2) = singlemat_x_mat(1:iend,j1,j2) & 244 244 & + A(j1,j3)*B(1:iend,j3,j2) … … 273 273 mat_x_singlemat = 0.0_jprb 274 274 275 doj2 = 1,m276 doj1 = 1,m277 doj3 = 1,m275 DO j2 = 1,m 276 DO j1 = 1,m 277 DO j3 = 1,m 278 278 mat_x_singlemat(1:iend,j1,j2) = mat_x_singlemat(1:iend,j1,j2) & 279 279 & + A(1:iend,j1,j3)*B(j3,j2) … … 311 311 312 312 identity_minus_mat_x_mat = - identity_minus_mat_x_mat 313 doj = 1,m313 DO j = 1,m 314 314 identity_minus_mat_x_mat(1:iend,j,j) & 315 315 & = 1.0_jprb + identity_minus_mat_x_mat(1:iend,j,j) … … 337 337 338 338 sparse_x_dense = 0.0_jprb 339 doj2 = 1,n2340 doj1 = 1,n1339 DO j2 = 1,n2 340 DO j1 = 1,n1 341 341 if (sparse(j1,j2) /= 0.0_jprb) then 342 342 sparse_x_dense(j1,:) = sparse_x_dense(j1,:) + sparse(j1,j2)*dense(j2,:) … … 376 376 mblock = m/3 377 377 m2block = 2*mblock 378 doj4 = 1,nrepeat378 DO j4 = 1,nrepeat 379 379 repeated_square = 0.0_jprb 380 380 ! Do the top-left (C, D, F & G) 381 doj2 = 1,m2block382 doj1 = 1,m2block383 doj3 = 1,m2block381 DO j2 = 1,m2block 382 DO j1 = 1,m2block 383 DO j3 = 1,m2block 384 384 repeated_square(j1,j2) = repeated_square(j1,j2) & 385 385 & + A(j1,j3)*A(j3,j2) … … 387 387 end do 388 388 end do 389 doj2 = m2block+1, m389 DO j2 = m2block+1, m 390 390 ! Do the top-right (E & H) 391 doj1 = 1,m2block392 doj3 = 1,m391 DO j1 = 1,m2block 392 DO j3 = 1,m 393 393 repeated_square(j1,j2) = repeated_square(j1,j2) & 394 394 & + A(j1,j3)*A(j3,j2) … … 396 396 end do 397 397 ! Do the bottom-right (I) 398 doj1 = m2block+1, m399 doj3 = m2block+1,m398 DO j1 = m2block+1, m 399 DO j3 = m2block+1,m 400 400 repeated_square(j1,j2) = repeated_square(j1,j2) & 401 401 & + A(j1,j3)*A(j3,j2) … … 409 409 else 410 410 ! Ordinary dense matrix 411 doj4 = 1,nrepeat411 DO j4 = 1,nrepeat 412 412 repeated_square = 0.0_jprb 413 doj2 = 1,m414 doj1 = 1,m415 doj3 = 1,m413 DO j2 = 1,m 414 DO j1 = 1,m 415 DO j3 = 1,m 416 416 repeated_square(j1,j2) = repeated_square(j1,j2) & 417 417 & + A(j1,j3)*A(j3,j2) … … 549 549 U33 = A(1:iend,3,3) - L31*A(1:iend,1,3) - L32*U23 550 550 551 doj = 1,3551 DO j = 1,3 552 552 ! Solve Ly = B(:,:,j) by forward substitution 553 553 ! y1 = B(:,1,j) … … 649 649 LU(1:iend,1:m,1:m) = A(1:iend,1:m,1:m) 650 650 651 doj2 = 1, m652 doj1 = 1, j2-1651 DO j2 = 1, m 652 DO j1 = 1, j2-1 653 653 s = LU(1:iend,j1,j2) 654 doj3 = 1, j1-1654 DO j3 = 1, j1-1 655 655 s = s - LU(1:iend,j1,j3) * LU(1:iend,j3,j2) 656 656 end do 657 657 LU(1:iend,j1,j2) = s 658 658 end do 659 doj1 = j2, m659 DO j1 = j2, m 660 660 s = LU(1:iend,j1,j2) 661 doj3 = 1, j2-1661 DO j3 = 1, j2-1 662 662 s = s - LU(1:iend,j1,j3) * LU(1:iend,j3,j2) 663 663 end do … … 666 666 if (j2 /= m) then 667 667 s = 1.0_jprb / LU(1:iend,j2,j2) 668 doj1 = j2+1, m668 DO j1 = j2+1, m 669 669 LU(1:iend,j1,j2) = LU(1:iend,j1,j2) * s 670 670 end do … … 691 691 692 692 ! First solve Ly=b 693 doj2 = 2, m694 doj1 = 1, j2-1693 DO j2 = 2, m 694 DO j1 = 1, j2-1 695 695 x(1:iend,j2) = x(1:iend,j2) - x(1:iend,j1)*LU(1:iend,j2,j1) 696 696 end do 697 697 end do 698 698 ! Now solve Ux=y 699 doj2 = m, 1, -1700 doj1 = j2+1, m699 DO j2 = m, 1, -1 700 DO j1 = j2+1, m 701 701 x(1:iend,j2) = x(1:iend,j2) - x(1:iend,j1)*LU(1:iend,j2,j1) 702 702 end do … … 722 722 call lu_factorization(n,iend,m,A,LU) 723 723 724 doj = 1, m724 DO j = 1, m 725 725 call lu_substitution(n,iend,m,LU,B(1:,1:m,j),X(1:iend,1:m,j)) 726 726 ! call lu_substitution(n,iend,m,LU,B(1:n,1:m,j),X(1:iend,1:m,j)) … … 835 835 836 836 ! Compute the 1-norms of A 837 doj3 = 1,m837 DO j3 = 1,m 838 838 sum_column(:) = 0.0_jprb 839 doj2 = 1,m840 doj1 = 1,iend839 DO j2 = 1,m 840 DO j1 = 1,iend 841 841 sum_column(j1) = sum_column(j1) + abs(A(j1,j2,j3)) 842 842 end do 843 843 end do 844 doj1 = 1,iend844 DO j1 = 1,iend 845 845 if (sum_column(j1) > normA(j1)) then 846 846 normA(j1) = sum_column(j1) … … 861 861 ! Scale the input matrices by a power of 2 862 862 scaling = 2.0_jprb**(-expo) 863 doj3 = 1,m864 doj2 = 1,m863 DO j3 = 1,m 864 DO j2 = 1,m 865 865 A(1:iend,j2,j3) = A(1:iend,j2,j3) * scaling 866 866 end do … … 872 872 873 873 V = c(8)*A6 + c(6)*A4 + c(4)*A2 874 doj3 = 1,m874 DO j3 = 1,m 875 875 V(:,j3,j3) = V(:,j3,j3) + c(2) 876 876 end do … … 878 878 V = c(7)*A6 + c(5)*A4 + c(3)*A2 879 879 ! Add a multiple of the identity matrix 880 doj3 = 1,m880 DO j3 = 1,m 881 881 V(:,j3,j3) = V(:,j3,j3) + c(1) 882 882 end do … … 887 887 888 888 ! Add the identity matrix 889 doj3 = 1,m889 DO j3 = 1,m 890 890 A(1:iend,j3,j3) = A(1:iend,j3,j3) + 1.0_jprb 891 891 end do 892 892 893 893 ! Loop through the matrices 894 doj1 = 1,iend894 DO j1 = 1,iend 895 895 if (expo(j1) > 0) then 896 896 ! Square matrix j1 expo(j1) times … … 1016 1016 1017 1017 ! Compute V * diag_rdivide_V 1018 doj1 = 1,31019 doj2 = 1,31018 DO j1 = 1,3 1019 DO j2 = 1,3 1020 1020 R(1:iend,j2,j1) = V(1:iend,j2,1)*diag_rdivide_V(1:iend,1,j1) & 1021 1021 & + 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 146 146 147 147 ! Loop through columns 148 dojcol = istartcol,iendcol148 DO jcol = istartcol,iendcol 149 149 150 150 ! Clear-sky calculation … … 203 203 is_clear_sky_layer = .true. 204 204 i_cloud_top = nlev+1 205 dojlev = 1,nlev205 DO jlev = 1,nlev 206 206 ! Compute combined gas+aerosol+cloud optical properties 207 207 if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then … … 212 212 end if 213 213 214 dojg = 1,ng214 DO jg = 1,ng 215 215 od_cloud_new(jg) = od_scaling(jg,jlev) & 216 216 & * od_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) … … 228 228 ! case that od_total > 0.0 and ssa_total > 0.0 but 229 229 ! od_total*ssa_total == 0 due to underflow 230 dojg = 1,ng230 DO jg = 1,ng 231 231 if (od_total(jg) > 0.0_jprb) then 232 232 scat_od_total(jg) = ssa(jg,jlev,jcol)*od(jg,jlev,jcol) & … … 247 247 else 248 248 249 dojg = 1,ng249 DO jg = 1,ng 250 250 if (od_total(jg) > 0.0_jprb) then 251 251 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 138 138 139 139 ! Loop through columns 140 dojcol = istartcol,iendcol140 DO jcol = istartcol,iendcol 141 141 ! Only perform calculation if sun above the horizon 142 142 if (single_level%cos_sza(jcol) > 0.0_jprb) then … … 155 155 else 156 156 ! Apply delta-Eddington scaling to the aerosol-gas mixture 157 dojlev = 1,nlev157 DO jlev = 1,nlev 158 158 od_total = od(:,jlev,jcol) 159 159 ssa_total = ssa(:,jlev,jcol) … … 205 205 if (total_cloud_cover >= config%cloud_fraction_threshold) then 206 206 ! Total-sky calculation 207 dojlev = 1,nlev207 DO jlev = 1,nlev 208 208 ! Compute combined gas+aerosol+cloud optical properties 209 209 if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then 210 dojg = 1,ng210 DO jg = 1,ng 211 211 od_cloud_new(jg) = od_scaling(jg,jlev) & 212 212 & * 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 164 164 integer :: jlev 165 165 166 dojlev = 1,nlev166 DO jlev = 1,nlev 167 167 ! The fraction of the total optical depth in the current layer 168 168 ! is proportional to the fraction of the mass of the atmosphere … … 192 192 & = StefanBoltzmann * single_level%skin_temperature(istartcol:iendcol)**4 & 193 193 & * single_level%lw_emissivity(istartcol:iendcol,1) 194 dojlev = 1,nlev+1194 DO jlev = 1,nlev+1 195 195 planck_hl(1,jlev,istartcol:iendcol) = StefanBoltzmann * thermodynamics%temperature_hl(istartcol:iendcol,jlev)**4 196 196 end do … … 200 200 & single_level%skin_temperature(istartcol:iendcol)) & 201 201 & * single_level%lw_emissivity(istartcol:iendcol,1) 202 dojlev = 1,nlev+1202 DO jlev = 1,nlev+1 203 203 planck_hl(1,jlev,istartcol:iendcol) = Pi*planck_function(config%mono_lw_wavelength, & 204 204 & thermodynamics%temperature_hl(istartcol:iendcol,jlev)) … … 261 261 ! Convert cloud mixing ratio into liquid and ice water path 262 262 ! in each layer 263 dojlev = 1, nlev264 dojcol = istartcol, iendcol263 DO jlev = 1, nlev 264 DO jcol = istartcol, iendcol 265 265 ! Factor to convert from gridbox-mean mass mixing ratio to 266 266 ! in-cloud water path involves the pressure difference in … … 288 288 289 289 if (config%iverbose >= 4) then 290 dojcol = istartcol,iendcol290 DO jcol = istartcol,iendcol 291 291 write(*,'(a,i0,a,f7.3,a,f7.3)') 'Profile ', jcol, ': shortwave optical depth = ', & 292 292 & sum(od_sw_cloud(1,:,jcol)*cloud%fraction(jcol,:)), & -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_overlap.F90
r4773 r5158 94 94 ! multiply by "factor" 95 95 denominator = 1.0_jprb 96 dojreg = 1,nreg96 DO jreg = 1,nreg 97 97 op_x_frac_min(jreg) = op(jreg) & 98 98 & * min(frac_upper(jreg), frac_lower(jreg)) … … 103 103 factor = 1.0_jprb / denominator 104 104 ! Create the random part of the overlap matrix 105 dojupper = 1,nreg106 dojlower = 1,nreg105 DO jupper = 1,nreg 106 DO jlower = 1,nreg 107 107 overlap_matrix(jupper,jlower) = factor & 108 108 & * (frac_lower(jlower)-op_x_frac_min(jlower)) & … … 115 115 116 116 ! Add on the maximum part of the overlap matrix 117 dojreg = 1,nreg117 DO jreg = 1,nreg 118 118 overlap_matrix(jreg,jreg) = overlap_matrix(jreg,jreg) & 119 119 & + op_x_frac_min(jreg) … … 368 368 369 369 ! Loop through each atmospheric column 370 dojcol = istartcol, iendcol370 DO jcol = istartcol, iendcol 371 371 ! For this column, outer space is treated as one clear-sky 372 372 ! region, so the fractions are assigned as such … … 381 381 ! half-level starting at 1 for the top-of-atmosphere, as well 382 382 ! as indexing each level starting at 1 for the top-most level. 383 dojlev = 1,nlev+1383 DO jlev = 1,nlev+1 384 384 ! Fraction of each region just below the interface 385 385 if (jlev > nlev) then … … 426 426 427 427 ! Convert to directional overlap matrices 428 dojupper = 1,nreg429 dojlower = 1,nreg428 DO jupper = 1,nreg 429 DO jlower = 1,nreg 430 430 if (frac_lower(jlower) >= frac_threshold) then 431 431 u_matrix(jupper,jlower,jlev,jcol) = overlap_matrix(jupper,jlower) & -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_pdf_sampler.F90
r4773 r5158 189 189 real(jprb) :: wfsd, wcdf 190 190 191 dojsamp = 1,nsamp191 DO jsamp = 1,nsamp 192 192 if (mask(jsamp)) then 193 193 ! Bilinear interpolation with bounds … … 237 237 real(jprb) :: wfsd, wcdf 238 238 239 dojz = 1,nz240 dojg = 1,ng239 DO jz = 1,nz 240 DO jg = 1,ng 241 241 if (cdf(jg, jz) > 0.0_jprb) then 242 242 ! Bilinear interpolation with bounds … … 291 291 real(jprb) :: wfsd, wcdf 292 292 293 dojz = 1,nz293 DO jz = 1,nz 294 294 295 295 if (mask(jz)) then 296 296 297 dojg = 1,ng297 DO jg = 1,ng 298 298 if (cdf(jg, jz) > 0.0_jprb) then 299 299 ! Bilinear interpolation with bounds -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_random_numbers.F90
r4773 r5158 165 165 ! populate the state 166 166 rseed = real(abs(this%iseed),jprd) 167 dojstr = 1,this%nmaxstreams167 DO jstr = 1,this%nmaxstreams 168 168 ! Note that nint returns an integer of type jpib (8-byte) 169 169 ! which may be converted to double if that is the type of … … 174 174 175 175 ! One warmup of the C++ minstd_rand algorithm 176 dojstr = 1,this%nmaxstreams176 DO jstr = 1,this%nmaxstreams 177 177 this%istate(jstr) = mod(IMinstdA * this%istate(jstr), IMinstdM) 178 178 end do … … 182 182 call random_seed(size=nseed) 183 183 allocate(iseednative(nseed)) 184 dojseed = 1,nseed184 DO jseed = 1,nseed 185 185 iseednative(jseed) = this%iseed + jseed - 1 186 186 end do … … 208 208 209 209 ! C++ minstd_rand algorithm 210 doi = 1, imax210 DO i = 1, imax 211 211 ! The following calculation is computed entirely with 8-byte 212 212 ! numbers (whether real or integer) … … 243 243 244 244 ! C++ minstd_ran algorithm 245 dojblock = 1,size(randnum,2)245 DO jblock = 1,size(randnum,2) 246 246 ! These lines should be vectorizable 247 doi = 1, imax247 DO i = 1, imax 248 248 this%istate(i) = mod(IMinstdA * this%istate(i), IMinstdM) 249 249 randnum(i,jblock) = IMinstdScale * this%istate(i) … … 278 278 279 279 ! C++ minstd_ran algorithm 280 dojblock = 1,size(randnum,2)280 DO jblock = 1,size(randnum,2) 281 281 if (mask(jblock)) then 282 282 ! These lines should be vectorizable 283 doi = 1, imax283 DO i = 1, imax 284 284 this%istate(i) = mod(IMinstdA * this%istate(i), IMinstdM) 285 285 randnum(i,jblock) = IMinstdScale * this%istate(i) … … 290 290 else 291 291 292 dojblock = 1,size(randnum,2)292 DO jblock = 1,size(randnum,2) 293 293 if (mask(jblock)) then 294 294 call random_number(randnum(:,jblock)) -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_regions.F90
r4773 r5158 122 122 ! sigma, then the 16th percentile of the lognormal is very 123 123 ! close to exp(mu-sigma). 124 dojcol = istartcol,iendcol125 dojlev = 1,nlev124 DO jcol = istartcol,iendcol 125 DO jlev = 1,nlev 126 126 if (cloud_fraction(jcol,jlev) < frac_threshold) then 127 127 reg_fracs(1,jlev,jcol) = 1.0_jprb … … 145 145 ! cover of the denser region - see region_fractions routine 146 146 ! below 147 dojcol = istartcol,iendcol148 dojlev = 1,nlev147 DO jcol = istartcol,iendcol 148 DO jlev = 1,nlev 149 149 if (cloud_fraction(jcol,jlev) < frac_threshold) then 150 150 reg_fracs(1,jlev,jcol) = 1.0_jprb -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_save.F90
r4853 r5158 1168 1168 & dim2_name="column", dim1_name="level", & 1169 1169 & units_str="1", long_name="Ozone mass mixing ratio") 1170 dojgas = 1,NMaxGases1170 DO jgas = 1,NMaxGases 1171 1171 if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then 1172 1172 write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr' … … 1259 1259 call gas%get(IO3, IMassMixingRatio, mixing_ratio) 1260 1260 call out_file%put("o3_mmr", mixing_ratio) 1261 dojgas = 1,NMaxGases1261 DO jgas = 1,NMaxGases 1262 1262 if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then 1263 1263 write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr' -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_single_level.F90
r4773 r5158 206 206 end if 207 207 208 dojcol = istartcol,iendcol208 DO jcol = istartcol,iendcol 209 209 this%iseed(jcol) = jcol 210 210 end do … … 275 275 276 276 sw_albedo_band = 0.0_jprb 277 dojband = 1,config%n_bands_sw278 dojalbedoband = 1,nalbedoband277 DO jband = 1,config%n_bands_sw 278 DO jalbedoband = 1,nalbedoband 279 279 if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then 280 dojcol = istartcol,iendcol280 DO jcol = istartcol,iendcol 281 281 sw_albedo_band(jcol,jband) & 282 282 & = sw_albedo_band(jcol,jband) & … … 292 292 if (allocated(this%sw_albedo_direct)) then 293 293 sw_albedo_band = 0.0_jprb 294 dojband = 1,config%n_bands_sw295 dojalbedoband = 1,nalbedoband294 DO jband = 1,config%n_bands_sw 295 DO jalbedoband = 1,nalbedoband 296 296 if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then 297 297 sw_albedo_band(istartcol:iendcol,jband) & … … 342 342 end if 343 343 lw_albedo_band = 0.0_jprb 344 dojband = 1,config%n_bands_lw345 dojalbedoband = 1,nalbedoband344 DO jband = 1,config%n_bands_lw 345 DO jalbedoband = 1,nalbedoband 346 346 if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then 347 dojcol = istartcol,iendcol347 DO jcol = istartcol,iendcol 348 348 lw_albedo_band(jcol,jband) & 349 349 & = lw_albedo_band(jcol,jband) & -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_lw.F90
r4773 r5158 313 313 314 314 ! Main loop over columns 315 dojcol = istartcol, iendcol315 DO jcol = istartcol, iendcol 316 316 ! -------------------------------------------------------- 317 317 ! Section 2: Prepare column-specific variables and arrays … … 325 325 ! cloud%crop_cloud_fraction has already been called 326 326 is_clear_sky_layer = .true. 327 dojlev = 1,nlev327 DO jlev = 1,nlev 328 328 if (cloud%fraction(jcol,jlev) > 0.0_jprb) then 329 329 is_clear_sky_layer(jlev) = .false. … … 336 336 ! In this section the reflectance, transmittance and sources 337 337 ! are computed for each layer 338 dojlev = 1,nlev ! Start at top-of-atmosphere338 DO jlev = 1,nlev ! Start at top-of-atmosphere 339 339 ! -------------------------------------------------------- 340 340 ! Section 3.1: Layer-specific initialization … … 395 395 ! approximate order of gas optical depth. 396 396 ng3D = ng 397 dojg = 1, ng397 DO jg = 1, ng 398 398 if (od_region(jg,1) > config%max_gas_od_3D) then 399 399 ng3D = jg-1 … … 488 488 end if 489 489 490 dojreg = 1, nreg-1490 DO jreg = 1, nreg-1 491 491 ! Compute lateral transfer rate from region jreg to 492 492 ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but … … 536 536 ! g-point, mapping from the cloud properties 537 537 ! defined in each band. 538 dojg = 1,ng538 DO jg = 1,ng 539 539 ! Mapping from g-point to band 540 540 iband = config%i_band_from_reordered_g_lw(jg) … … 544 544 545 545 ! Loop over cloudy regions 546 dojreg = 2,nreg546 DO jreg = 2,nreg 547 547 ! Add scaled cloud optical depth to clear-sky value 548 548 od_region(jg,jreg) = od_region(jg,1) & … … 603 603 ! exponential, then computing the various 604 604 ! transmission/reflectance matrices from that. 605 dojreg = 1,nregActive605 DO jreg = 1,nregActive 606 606 ! Write the diagonal elements of -Gamma1*z1 607 607 Gamma_z1(1:ng3D,jreg,jreg) & … … 628 628 ! non-zeros in the cloudy-sky regions even if the cloud 629 629 ! fraction is zero 630 dojreg = nregActive+1,nreg630 DO jreg = nregActive+1,nreg 631 631 Gamma_z1(1:ng3D,jreg,jreg) = Gamma_z1(1:ng3D,1,1) 632 632 Gamma_z1(1:ng3D,nreg+jreg,jreg) = Gamma_z1(1:ng3D,nreg+1,1) … … 656 656 end if 657 657 658 dojreg = 1,nregActive-1658 DO jreg = 1,nregActive-1 659 659 ! Add the terms assocated with 3D transport to Gamma1*z1. 660 660 ! First the rate from region jreg to jreg+1 … … 775 775 ! formulas for reflectance and transmittance, and equivalent 776 776 ! solutions to the coupled ODEs for the sources. 777 dojreg = 2, nregActive777 DO jreg = 2, nregActive 778 778 call calc_reflectance_transmittance_lw(ng-ng3D, & 779 779 & od_region(ng3D+1:ng,jreg), & … … 806 806 ! Calculate the upwelling radiation emitted by the surface, and 807 807 ! copy the surface albedo into total_albedo 808 dojreg = 1,nreg809 dojg = 1,ng808 DO jreg = 1,nreg 809 DO jg = 1,ng 810 810 ! region_fracs(jreg,nlev,jcol) is the fraction of each 811 811 ! region in the lowest model level … … 816 816 ! Equivalent surface values for computing clear-sky fluxes 817 817 if (config%do_clear) then 818 dojg = 1,ng818 DO jg = 1,ng 819 819 total_source_clear(jg,nlev+1) = emission(jg,jcol) 820 820 end do … … 826 826 ! Loop back up through the atmosphere computing the total albedo 827 827 ! and the total upwelling due to emission below each level 828 dojlev = nlev,1,-1828 DO jlev = nlev,1,-1 829 829 if (config%do_clear) then 830 830 ! Use adding method for clear-sky arrays; note that there … … 876 876 total_albedo_below = 0.0_jprb 877 877 total_source_below = 0.0_jprb 878 dojreg = 1,nreg878 DO jreg = 1,nreg 879 879 inv_denom_scalar(:) = 1.0_jprb / (1.0_jprb & 880 880 & - total_albedo(:,jreg,jreg,jlev+1)*reflectance(:,jreg,jreg,jlev)) … … 918 918 ! essentially diag(total_albedo) = matmul(transpose(v_matrix), 919 919 ! diag(total_albedo_below)). 920 dojreg = 1,nreg921 dojreg2 = 1,nreg920 DO jreg = 1,nreg 921 DO jreg2 = 1,nreg 922 922 total_albedo(:,jreg,jreg,jlev) & 923 923 & = total_albedo(:,jreg,jreg,jlev) & … … 966 966 967 967 ! Final loop back down through the atmosphere to compute fluxes 968 dojlev = 1,nlev968 DO jlev = 1,nlev 969 969 if (config%do_clear) then 970 970 ! Scalar operations for clear-sky fluxes … … 999 999 & flux_dn_above) + total_source(:,:,jlev+1) 1000 1000 else 1001 dojreg = 1,nreg1001 DO jreg = 1,nreg 1002 1002 ! Scalar operations for all regions, requiring that 1003 1003 ! reflectance, transmittance and total_albedo are diagonal -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spartacus_sw.F90
r4853 r5158 328 328 329 329 ! Main loop over columns 330 dojcol = istartcol, iendcol330 DO jcol = istartcol, iendcol 331 331 ! -------------------------------------------------------- 332 332 ! Section 2: Prepare column-specific variables and arrays … … 404 404 is_clear_sky_layer = .true. 405 405 i_cloud_top = nlev+1 406 dojlev = nlev,1,-1406 DO jlev = nlev,1,-1 407 407 if (cloud%fraction(jcol,jlev) > 0.0_jprb) then 408 408 is_clear_sky_layer(jlev) = .false. … … 416 416 ! In this section the reflectance, transmittance and sources 417 417 ! are computed for each layer 418 dojlev = 1,nlev ! Start at top-of-atmosphere418 DO jlev = 1,nlev ! Start at top-of-atmosphere 419 419 ! -------------------------------------------------------- 420 420 ! Section 3.1: Layer-specific initialization … … 469 469 ! approximate order of gas optical depth. 470 470 ng3D = ng 471 dojg = 1, ng471 DO jg = 1, ng 472 472 if (od_region(jg,1) > config%max_gas_od_3D) then 473 473 ng3D = jg-1 … … 553 553 end if 554 554 555 dojreg = 1, nreg-1555 DO jreg = 1, nreg-1 556 556 ! Compute lateral transfer rates from region jreg to 557 557 ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but … … 615 615 ! g-point, mapping from the cloud properties 616 616 ! defined in each band. 617 dojg = 1,ng617 DO jg = 1,ng 618 618 ! Mapping from g-point to band 619 619 iband = config%i_band_from_reordered_g_sw(jg) … … 629 629 630 630 ! Loop over cloudy regions 631 dojreg = 2,nreg631 DO jreg = 2,nreg 632 632 scat_od_cloud = od_cloud(iband,jlev,jcol) & 633 633 & * ssa_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol) … … 679 679 ! exponential, then computing the various 680 680 ! transmission/reflectance matrices from that. 681 dojreg = 1,nregactive681 DO jreg = 1,nregactive 682 682 ! Write the diagonal elements of -Gamma1*z1 683 683 Gamma_z1(1:ng3D,jreg,jreg) & … … 701 701 end do 702 702 703 dojreg = 1,nregactive-1703 DO jreg = 1,nregactive-1 704 704 ! Write the elements of -Gamma1*z1 concerned with 3D 705 705 ! transport … … 820 820 ! Compute reflectance, transmittance and associated terms 821 821 ! for each cloudy region, using the Meador-Weaver formulas 822 dojreg = 2, nregactive822 DO jreg = 2, nregactive 823 823 call calc_reflectance_transmittance_sw(ng-ng3D, & 824 824 & mu0, & … … 853 853 ! beam incident on the surface, and copy the surface albedo 854 854 ! into total_albedo 855 dojreg = 1,nreg856 dojg = 1,ng855 DO jreg = 1,nreg 856 DO jg = 1,ng 857 857 total_albedo(jg,jreg,jreg,nlev+1) = albedo_diffuse(jg,jcol) 858 858 total_albedo_direct(jg,jreg,jreg,nlev+1) & … … 875 875 ! Work back up through the atmosphere computing the total albedo 876 876 ! of the atmosphere below that point using the adding method 877 dojlev = nlev,1,-1877 DO jlev = nlev,1,-1 878 878 879 879 ! -------------------------------------------------------- … … 995 995 ! radiation: 996 996 total_albedo(:,:,:,jlev) = 0.0_jprb 997 dojreg = 1,nreg ! Target layer (jlev-1)998 dojreg2 = 1,nreg ! Current layer (jlev)997 DO jreg = 1,nreg ! Target layer (jlev-1) 998 DO jreg2 = 1,nreg ! Current layer (jlev) 999 999 total_albedo(:,jreg,jreg,jlev) = total_albedo(:,jreg,jreg,jlev) & 1000 1000 & + sum(total_albedo_below(:,:,jreg2),2) & … … 1004 1004 ! ...then direct radiation: 1005 1005 total_albedo_direct(:,:,:,jlev) = 0.0_jprb 1006 dojreg = 1,nreg ! Target layer (jlev-1)1007 dojreg2 = 1,nreg ! Current layer (jlev)1006 DO jreg = 1,nreg ! Target layer (jlev-1) 1007 DO jreg2 = 1,nreg ! Current layer (jlev) 1008 1008 total_albedo_direct(:,jreg,jreg,jlev) = total_albedo_direct(:,jreg,jreg,jlev) & 1009 1009 & + sum(total_albedo_below_direct(:,:,jreg2),2) & … … 1031 1031 ! First diffuse radiation: 1032 1032 albedo_part = total_albedo_below 1033 dojreg = 1,nreg1033 DO jreg = 1,nreg 1034 1034 albedo_part(:,jreg,jreg) = 0.0_jprb 1035 1035 end do … … 1040 1040 ! ...then direct radiation: 1041 1041 albedo_part = total_albedo_below_direct 1042 dojreg = 1,nreg1042 DO jreg = 1,nreg 1043 1043 albedo_part(:,jreg,jreg) = 0.0_jprb 1044 1044 end do … … 1059 1059 ! essentially diag(total_albedo) += matmul(transpose(v_matrix), 1060 1060 ! diag(total_albedo_below)). 1061 dojreg = 1,nreg1062 dojreg2 = 1,nreg1061 DO jreg = 1,nreg 1062 DO jreg2 = 1,nreg 1063 1063 total_albedo(:,jreg,jreg,jlev) & 1064 1064 & = total_albedo(:,jreg,jreg,jlev) & … … 1075 1075 ! "Explicit entrapment" 1076 1076 1077 dojreg2 = 1,nreg1077 DO jreg2 = 1,nreg 1078 1078 ! Loop through each region in the lower layer. For one 1079 1079 ! of the regions in the lower layer, we are imagining it … … 1110 1110 & / max(config%cloud_fraction_threshold, region_fracs(jreg2,jlev,jcol)) 1111 1111 1112 dojreg = 1, nreg-11112 DO jreg = 1, nreg-1 1113 1113 ! Compute lateral transfer rates from region jreg to 1114 1114 ! jreg+1 as before, but without length scale which … … 1139 1139 inv_effective_size = min(cloud%inv_cloud_effective_size(jcol,jlev-1), & 1140 1140 & 1.0_jprb/config%min_cloud_effective_size) 1141 dojreg = 1,nreg-11141 DO jreg = 1,nreg-1 1142 1142 ! Diffuse transport down and up with one random 1143 1143 ! scattering event … … 1167 1167 ! trival limit, so we cap the maximum input to expm by 1168 1168 ! scaling down if necessary 1169 dojg = 1,ng1169 DO jg = 1,ng 1170 1170 max_entr = -min(entrapment(jg,1,1),entrapment(jg,2,2)) 1171 1171 if (max_entr > config%max_cloud_od) then … … 1197 1197 #ifndef EXPLICIT_EDGE_ENTRAPMENT 1198 1198 ! Scale to get the contribution to the diffuse albedo 1199 dojreg3 = 1,nreg1200 dojreg = 1,nreg1199 DO jreg3 = 1,nreg 1200 DO jreg = 1,nreg 1201 1201 albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) & 1202 1202 & * v_matrix(jreg2,jreg,jlev,jcol) * total_albedo_below(:,jreg2,jreg2) … … 1210 1210 albedo_part = 0.0_jprb 1211 1211 ! Scale to get the contribution to the diffuse albedo 1212 dojreg3 = 1,nreg ! TO upper region1213 dojreg = 1,nreg ! FROM upper region1212 DO jreg3 = 1,nreg ! TO upper region 1213 DO jreg = 1,nreg ! FROM upper region 1214 1214 transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) & 1215 1215 & * cloud%overlap_param(jcol,jlev-1) & 1216 1216 & * min(region_fracs(jreg,jlev,jcol), region_fracs(jreg,jlev,jcol)) & 1217 1217 & / max(config%cloud_fraction_threshold, region_fracs(jreg,jlev,jcol)) 1218 dojreg4 = 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) 1219 1219 if (.not. (jreg4 == jreg .and. jreg4 /= jreg2)) then 1220 1220 albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) + entrapment(:,jreg3,jreg) & … … 1236 1236 ! Now do the same for the direct albedo 1237 1237 entrapment = 0.0_jprb 1238 dojreg = 1,nreg-11238 DO jreg = 1,nreg-1 1239 1239 ! Direct transport down and diffuse up with one random 1240 1240 ! scattering event … … 1264 1264 ! trival limit, so we cap the maximum input to expm by 1265 1265 ! scaling down if necessary 1266 dojg = 1,ng1266 DO jg = 1,ng 1267 1267 max_entr = -min(entrapment(jg,1,1),entrapment(jg,2,2)) 1268 1268 if (max_entr > config%max_cloud_od) then … … 1289 1289 1290 1290 #ifndef EXPLICIT_EDGE_ENTRAPMENT 1291 dojreg3 = 1,nreg1292 dojreg = 1,nreg1291 DO jreg3 = 1,nreg 1292 DO jreg = 1,nreg 1293 1293 albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) & 1294 1294 & * v_matrix(jreg2,jreg,jlev,jcol) * total_albedo_below_direct(:,jreg2,jreg2) … … 1298 1298 entrapment = albedo_part 1299 1299 albedo_part = 0.0_jprb 1300 dojreg3 = 1,nreg1301 dojreg = 1,nreg1300 DO jreg3 = 1,nreg 1301 DO jreg = 1,nreg 1302 1302 transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) & 1303 1303 & * cloud%overlap_param(jcol,jlev-1) & 1304 1304 & * min(region_fracs(jreg,jlev,jcol), region_fracs(jreg,jlev-1,jcol)) & 1305 1305 & / max(config%cloud_fraction_threshold, region_fracs(jreg,jlev,jcol)) 1306 dojreg4 = 1,nreg1306 DO jreg4 = 1,nreg 1307 1307 if (.not. (jreg4 == jreg .and. jreg4 /= jreg2)) then 1308 1308 albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) + entrapment(:,jreg3,jreg) & … … 1343 1343 end if 1344 1344 1345 dojreg = 1,nreg ! Target layer (jlev-1)1346 dojreg2 = 1,nregactive ! Current layer (jlev)1345 DO jreg = 1,nreg ! Target layer (jlev-1) 1346 DO jreg2 = 1,nregactive ! Current layer (jlev) 1347 1347 x_direct_above(:,jreg) = x_direct_above(:,jreg) & 1348 1348 & + x_direct(:,jreg2) * v_matrix(jreg2,jreg,jlev,jcol) … … 1368 1368 ! Direct downwelling flux (into a plane perpendicular to the 1369 1369 ! sun) entering the top of each region in the top-most layer 1370 dojreg = 1,nreg1370 DO jreg = 1,nreg 1371 1371 direct_dn_below(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol) 1372 1372 end do … … 1421 1421 1422 1422 ! Final loop back down through the atmosphere to compute fluxes 1423 dojlev = 1,nlev1423 DO jlev = 1,nlev 1424 1424 1425 1425 #ifdef PRINT_ENTRAPMENT_DATA … … 1675 1675 & + tan_diffuse_angle_3d*tan_diffuse_angle_3d) * 0.5_jprb 1676 1676 1677 dojreg = istartreg,iendreg1677 DO jreg = istartreg,iendreg 1678 1678 ! Geometric series enhancement due to multiple scattering: the 1679 1679 ! amplitude enhancement is equal to the limit of -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_spectral_definition.F90
r4853 r5158 201 201 else 202 202 find_wavenumber = 1 203 dowhile (wavenumber > this%wavenumber2(find_wavenumber) &203 DO while (wavenumber > this%wavenumber2(find_wavenumber) & 204 204 & .and. find_wavenumber < this%nwav) 205 205 find_wavenumber = find_wavenumber + 1 … … 284 284 end if 285 285 286 dojband = 1,this%nband286 DO jband = 1,this%nband 287 287 weight = 0.0_jprb 288 dojwav = 1,nwav288 DO jwav = 1,nwav 289 289 ! Work out wavenumber range for which this cloud wavenumber 290 290 ! will be applicable … … 321 321 ! Find interpolating points 322 322 iwav = 2 323 dowhile (wavenumber(iwav) < this%wavenumber2_band(jband))323 DO while (wavenumber(iwav) < this%wavenumber2_band(jband)) 324 324 iwav = iwav+1 325 325 end do … … 358 358 mapping = 0.0_jprb 359 359 ! Loop over wavenumbers representing cloud 360 dojwav = 1,nwav360 DO jwav = 1,nwav 361 361 ! Clear the weights. The weight says for one wavenumber in the 362 362 ! cloud file, what is its fractional contribution to each of … … 400 400 & / (this%wavenumber2(isd1)-this%wavenumber1(isd1)) 401 401 if (isd1-isd0 > 1) then 402 doisd = isd0+1,isd1-1402 DO isd = isd0+1,isd1-1 403 403 ! Intermediate trapezia 404 404 weight(isd) = 0.5_jprb * (this%wavenumber1(isd)+this%wavenumber2(isd) & … … 447 447 & / (this%wavenumber2(isd1)-this%wavenumber1(isd1)) 448 448 if (isd2-isd1 > 1) then 449 doisd = isd1+1,isd2-1449 DO isd = isd1+1,isd2-1 450 450 ! Intermediate trapezia 451 451 weight(isd) = weight(isd) + 0.5_jprb * (2.0_jprb*wavenum2 & … … 468 468 weight = weight * planck_weight 469 469 470 dojg = 1,this%ng470 DO jg = 1,this%ng 471 471 mapping(jg, jwav) = sum(weight * this%gpoint_fraction(:,jg)) 472 472 end do … … 478 478 479 479 ! Normalize mapping matrix 480 dojg = 1,this%ng480 DO jg = 1,this%ng 481 481 mapping(jg,:) = mapping(jg,:) * (1.0_jprb/sum(mapping(jg,:))) 482 482 end do … … 589 589 ! Check wavelength is monotonically increasing 590 590 if (ninterval > 2) then 591 dojint = 2,ninterval-1591 DO jint = 2,ninterval-1 592 592 if (wavelength_bound(jint) <= wavelength_bound(jint-1)) then 593 593 write(nulerr, '(a)') '*** Error: wavelength bounds must be monotonically increasing' … … 609 609 end if 610 610 611 dojband = 1,this%nband612 dojint = 1,ninterval611 DO jband = 1,this%nband 612 DO jint = 1,ninterval 613 613 if (jint == 1) then 614 614 ! First input interval in wavelength space: lower … … 692 692 693 693 ! All bounded intervals 694 dojint = 2,ninterval-1694 DO jint = 2,ninterval-1 695 695 wavenumber1_bound = 0.01_jprb / wavelength_bound(jint) 696 696 wavenumber2_bound = 0.01_jprb / wavelength_bound(jint-1) … … 710 710 end if 711 711 712 dojg = 1,this%ng713 dojin = 1,ninput712 DO jg = 1,this%ng 713 DO jin = 1,ninput 714 714 mapping(jin,jg) = sum(this%gpoint_fraction(:,jg) * planck, & 715 715 & mask=(i_input==jin)) … … 723 723 724 724 ! Loop through all intervals 725 dojint = 1,ninterval725 DO jint = 1,ninterval 726 726 ! Loop through the wavenumbers for gpoint_fraction 727 dojwav = 1,this%nwav727 DO jwav = 1,this%nwav 728 728 if (jint == 1) then 729 729 ! First input interval in wavelength space: lower … … 758 758 end do 759 759 if (use_fluxes_local) then 760 dojg = 1,this%ng760 DO jg = 1,this%ng 761 761 mapping(:,jg) = mapping(:,jg) / sum(this%gpoint_fraction(:,jg) * planck) 762 762 end do … … 769 769 if (.not. use_fluxes_local) then 770 770 ! Normalize mapping matrix 771 dojg = 1,size(mapping,dim=2)771 DO jg = 1,size(mapping,dim=2) 772 772 mapping(:,jg) = mapping(:,jg) * (1.0_jprb/sum(mapping(:,jg))) 773 773 end do … … 828 828 is_band_unassigned = .true. 829 829 830 dojint = 1,ninterval830 DO jint = 1,ninterval 831 831 inext = minloc(wavelength1, dim=1, mask=is_band_unassigned) 832 832 if (jint > 1) then … … 873 873 write(nulout, '(a,i0,a,i0,a)') ' Mapping from ', nin, ' values to ', nout, ' bands (wavenumber ranges in cm-1)' 874 874 if (nout <= 40) then 875 dojout = 1,nout875 DO jout = 1,nout 876 876 write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(jout)), ' to', & 877 877 & nint(this%wavenumber2_band(jout)), ':' 878 dojin = 1,nin878 DO jin = 1,nin 879 879 write(nulout,'(f5.2)',advance='no') mapping(jin,jout) 880 880 end do … … 882 882 end do 883 883 else 884 dojout = 1,30884 DO jout = 1,30 885 885 write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(jout)), ' to', & 886 886 & nint(this%wavenumber2_band(jout)), ':' 887 dojin = 1,nin887 DO jin = 1,nin 888 888 write(nulout,'(f5.2)',advance='no') mapping(jin,jout) 889 889 end do … … 893 893 write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(nout)), ' to', & 894 894 & nint(this%wavenumber2_band(nout)), ':' 895 dojin = 1,nin895 DO jin = 1,nin 896 896 write(nulout,'(f5.2)',advance='no') mapping(jin,nout) 897 897 end do … … 901 901 write(nulout, '(a,i0,a,i0,a)') ' Mapping from ', nin, ' values to ', nout, ' g-points' 902 902 if (nout <= 40) then 903 dojout = 1,nout903 DO jout = 1,nout 904 904 write(nulout,'(i3,a)',advance='no') jout, ':' 905 dojin = 1,nin905 DO jin = 1,nin 906 906 write(nulout,'(f5.2)',advance='no') mapping(jin,jout) 907 907 end do … … 909 909 end do 910 910 else 911 dojout = 1,30911 DO jout = 1,30 912 912 write(nulout,'(i3,a)',advance='no') jout, ':' 913 dojin = 1,nin913 DO jin = 1,nin 914 914 write(nulout,'(f5.2)',advance='no') mapping(jin,jout) 915 915 end do … … 918 918 write(nulout,'(a)') ' ...' 919 919 write(nulout,'(i3,a)',advance='no') nout, ':' 920 dojin = 1,nin920 DO jin = 1,nin 921 921 write(nulout,'(f5.2)',advance='no') mapping(jin,nout) 922 922 end do -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/radiation/radiation_thermodynamics.F90
r4773 r5158 143 143 end if 144 144 145 dojlev = 1,nlev146 dojcol = istartcol,iendcol145 DO jlev = 1,nlev 146 DO jcol = istartcol,iendcol 147 147 pressure = 0.5 * (this%pressure_hl(jcol,jlev)+this%pressure_hl(jcol,jlev+1)) 148 148 temperature = 0.5 * (this%temperature_hl(jcol,jlev)+this%temperature_hl(jcol,jlev+1)) … … 265 265 ! be half the separation between the half-levels at the edge of 266 266 ! the two adjacent layers 267 dojlev = 2,nlev-1267 DO jlev = 2,nlev-1 268 268 layer_separation(i1:i2,jlev) = (0.5_jprb * R_over_g) * temperature_hl(i1:i2,jlev+1) & 269 269 & * log(pressure_hl(i1:i2,jlev+2)/pressure_hl(i1:i2,jlev)) … … 276 276 ! don't take the logarithm of the last pressure in each column. 277 277 278 dojlev = 1,nlev-2278 DO jlev = 1,nlev-2 279 279 layer_separation(i1:i2,jlev) = (0.5_jprb * R_over_g) * temperature_hl(i1:i2,jlev+1) & 280 280 & * 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 202 202 203 203 ! Main loop over columns 204 dojcol = istartcol, iendcol204 DO jcol = istartcol, iendcol 205 205 ! -------------------------------------------------------- 206 206 ! Section 2: Prepare column-specific variables and arrays … … 211 211 is_clear_sky_layer = .true. 212 212 i_cloud_top = nlev+1 213 dojlev = 1,nlev213 DO jlev = 1,nlev 214 214 if (cloud%fraction(jcol,jlev) > 0.0_jprb) then 215 215 is_clear_sky_layer(jlev) = .false. … … 268 268 ! Save the spectral fluxes if required 269 269 if (config%do_save_spectral_flux) then 270 dojlev = 1,nlev+1270 DO jlev = 1,nlev+1 271 271 call indexed_sum(flux_up_clear(:,jlev), & 272 272 & config%i_spec_from_reordered_g_lw, & … … 294 294 end if 295 295 296 dojlev = i_cloud_top,nlev ! Start at cloud top and work down296 DO jlev = i_cloud_top,nlev ! Start at cloud top and work down 297 297 298 298 ! Copy over clear-sky properties … … 308 308 source_dn(:,2:,jlev) = 0.0_jprb 309 309 else 310 dojreg = 2,nreg310 DO jreg = 2,nreg 311 311 ! Cloudy sky 312 312 ! Add scaled cloud optical depth to clear-sky value … … 358 358 end do 359 359 ! Emission is scaled by the size of each region 360 dojreg = 1,nregions360 DO jreg = 1,nregions 361 361 source_up(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * source_up(:,jreg,jlev) 362 362 source_dn(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * source_dn(:,jreg,jlev) … … 375 375 ! Calculate the upwelling radiation emitted by the surface, and 376 376 ! copy the surface albedo into total_albedo 377 dojreg = 1,nregions378 dojg = 1,ng377 DO jreg = 1,nregions 378 DO jg = 1,ng 379 379 ! region_fracs(jreg,nlev,jcol) is the fraction of each region in the 380 380 ! lowest model level … … 387 387 ! atmosphere and the total upwelling due to emission below each 388 388 ! level below using the adding method 389 dojlev = nlev,i_cloud_top,-1389 DO jlev = nlev,i_cloud_top,-1 390 390 391 391 total_albedo_below = 0.0_jprb … … 394 394 total_albedo_below = 0.0_jprb 395 395 total_source_below = 0.0_jprb 396 dojg = 1,ng396 DO jg = 1,ng 397 397 inv_denom(jg,1) = 1.0_jprb & 398 398 & / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance(jg,1,jlev)) … … 428 428 ! the operation we perform is essentially diag(total_albedo) 429 429 ! = matmul(transpose(v_matrix), diag(total_albedo_below)). 430 dojreg = 1,nregions431 dojreg2 = 1,nregions430 DO jreg = 1,nregions 431 DO jreg2 = 1,nregions 432 432 total_albedo(:,jreg,jlev) & 433 433 & = total_albedo(:,jreg,jlev) & … … 445 445 ! Section 6: Copy over downwelling fluxes above cloud top 446 446 ! -------------------------------------------------------- 447 dojlev = 1,i_cloud_top447 DO jlev = 1,i_cloud_top 448 448 if (config%do_clear) then 449 449 ! Clear-sky fluxes have already been averaged: use these … … 476 476 & flux%lw_up_band(:,jcol,i_cloud_top)) 477 477 end if 478 dojlev = i_cloud_top-1,1,-1478 DO jlev = i_cloud_top-1,1,-1 479 479 flux_up(:,1) = trans_clear(:,jlev)*flux_up(:,1) + source_up_clear(:,jlev) 480 480 flux%lw_up(jcol,jlev) = sum(flux_up(:,1)) … … 494 494 ! scattering layer, using overlap matrix to translate to the 495 495 ! regions of the first layer of cloud 496 dojreg = 1,nregions496 DO jreg = 1,nregions 497 497 flux_dn(:,jreg) = v_matrix(jreg,1,i_cloud_top,jcol)*flux_dn_clear(:,i_cloud_top) 498 498 end do 499 499 500 500 ! Final loop back down through the atmosphere to compute fluxes 501 dojlev = i_cloud_top,nlev501 DO jlev = i_cloud_top,nlev 502 502 503 503 if (is_clear_sky_layer(jlev)) then 504 dojg = 1,ng504 DO jg = 1,ng 505 505 flux_dn(jg,1) = (transmittance(jg,1,jlev)*flux_dn(jg,1) & 506 506 & + 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 127 127 ! entire clear-sky atmosphere at once 128 128 real(jprb), dimension(config%n_g_sw, nlev) & 129 129 & :: reflectance_clear, transmittance_clear, & 130 130 & ref_dir_clear, trans_dir_diff_clear, trans_dir_dir_clear 131 131 … … 198 198 199 199 ! Main loop over columns 200 dojcol = istartcol, iendcol200 DO jcol = istartcol, iendcol 201 201 ! -------------------------------------------------------- 202 202 ! Section 2: Prepare column-specific variables and arrays … … 251 251 ! cloud%crop_cloud_fraction has already been called 252 252 is_clear_sky_layer = .true. 253 dojlev = 1,nlev253 DO jlev = 1,nlev 254 254 if (cloud%fraction(jcol,jlev) > 0.0_jprb) then 255 255 is_clear_sky_layer(jlev) = .false. … … 271 271 272 272 ! Cloudy layers 273 dojlev = 1,nlev ! Start at top-of-atmosphere273 DO jlev = 1,nlev ! Start at top-of-atmosphere 274 274 if (.not. is_clear_sky_layer(jlev)) then 275 dojreg = 2,nregions276 dojg = 1,ng275 DO jreg = 2,nregions 276 DO jg = 1,ng 277 277 ! Mapping from g-point to band 278 278 iband = config%i_band_from_reordered_g_sw(jg) … … 315 315 316 316 ! Copy surface albedos in clear-sky region 317 dojg = 1,ng317 DO jg = 1,ng 318 318 total_albedo(jg,1,nlev+1) = albedo_diffuse(jg,jcol) 319 319 total_albedo_direct(jg,1,nlev+1) & … … 324 324 ! underneath 325 325 if (.not. is_clear_sky_layer(nlev)) then 326 dojreg = 2,nregions326 DO jreg = 2,nregions 327 327 total_albedo(:,jreg,nlev+1) = total_albedo(:,1,nlev+1) 328 328 total_albedo_direct(:,jreg,nlev+1) = total_albedo_direct(:,1,nlev+1) … … 337 337 ! Work up from the surface computing the total albedo of the 338 338 ! atmosphere below that point using the adding method 339 dojlev = nlev,1,-1339 DO jlev = nlev,1,-1 340 340 341 341 total_albedo_below = 0.0_jprb … … 346 346 ! "below" quantities since with no cloud overlap to worry 347 347 ! about, these are the same 348 dojg = 1,ng348 DO jg = 1,ng 349 349 inv_denom(jg,1) = 1.0_jprb & 350 350 & / (1.0_jprb - total_albedo_clear(jg,jlev+1)*reflectance_clear(jg,jlev)) … … 361 361 362 362 ! All-sky fluxes: first the clear region 363 dojg = 1,ng363 DO jg = 1,ng 364 364 inv_denom(jg,1) = 1.0_jprb & 365 365 & / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance_clear(jg,jlev)) … … 375 375 ! Then the cloudy regions if any cloud is present in this layer 376 376 if (.not. is_clear_sky_layer(jlev)) then 377 dojreg = 2,nregions378 dojg = 1,ng377 DO jreg = 2,nregions 378 DO jg = 1,ng 379 379 inv_denom(jg,jreg) = 1.0_jprb / (1.0_jprb & 380 380 & - total_albedo(jg,jreg,jlev+1)*reflectance(jg,jreg,jlev)) … … 400 400 ! the operation we perform is essentially diag(total_albedo) 401 401 ! = matmul(transpose(v_matrix)), diag(total_albedo_below)). 402 dojreg = 1,nregions403 dojreg2 = 1,nregions402 DO jreg = 1,nregions 403 DO jreg2 = 1,nregions 404 404 total_albedo(:,jreg,jlev) & 405 405 & = total_albedo(:,jreg,jlev) & … … 426 426 ! Direct downwelling flux (into a plane perpendicular to the 427 427 ! sun) entering the top of each region in the top-most layer 428 dojreg = 1,nregions428 DO jreg = 1,nregions 429 429 direct_dn(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol) 430 430 end do … … 495 495 496 496 ! Final loop back down through the atmosphere to compute fluxes 497 dojlev = 1,nlev497 DO jlev = 1,nlev 498 498 if (config%do_clear) then 499 dojg = 1,ng499 DO jg = 1,ng 500 500 flux_dn_clear(jg) = (transmittance_clear(jg,jlev)*flux_dn_clear(jg) + direct_dn_clear(jg) & 501 501 & * (trans_dir_dir_clear(jg,jlev)*total_albedo_clear_direct(jg,jlev+1)*reflectance_clear(jg,jlev) & … … 509 509 510 510 ! All-sky fluxes: first the clear region... 511 dojg = 1,ng511 DO jg = 1,ng 512 512 flux_dn(jg,1) = (transmittance_clear(jg,jlev)*flux_dn(jg,1) + direct_dn(jg,1) & 513 513 & * (trans_dir_dir_clear(jg,jlev)*total_albedo_direct(jg,1,jlev+1)*reflectance_clear(jg,jlev) & … … 525 525 direct_dn(:,2:nregions)= 0.0_jprb 526 526 else 527 dojreg = 2,nregions528 dojg = 1,ng527 DO jreg = 2,nregions 528 DO jg = 1,ng 529 529 flux_dn(jg,jreg) = (transmittance(jg,jreg,jlev)*flux_dn(jg,jreg) + direct_dn(jg,jreg) & 530 530 & * (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 72 72 ! Added for DWD (2020) 73 73 !NEC$ shortloop 74 dojg = 1, ng74 DO jg = 1, ng 75 75 ! Fu et al. (1997), Eq 2.9 and 2.10: 76 76 ! gamma1(jg) = LwDiffusivity * (1.0_jprb - 0.5_jprb*ssa(jg) & … … 122 122 ! Added for DWD (2020) 123 123 !NEC$ shortloop 124 dojg = 1, ng124 DO jg = 1, ng 125 125 ! gamma1(jg) = 2.0_jprb - ssa(jg) * (1.25_jprb + 0.75_jprb*g(jg)) 126 126 ! gamma2(jg) = 0.75_jprb *(ssa(jg) * (1.0_jprb - g(jg))) … … 195 195 ! Added for DWD (2020) 196 196 !NEC$ shortloop 197 dojg = 1, ng197 DO jg = 1, ng 198 198 if (od(jg) > 1.0e-3_jprd) then 199 199 k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), & … … 291 291 #endif 292 292 293 dojg = 1, ng293 DO jg = 1, ng 294 294 factor = (LwDiffusivityWP * 0.5_jprb) * ssa(jg) 295 295 gamma1 = LwDiffusivityWP - factor*(1.0_jprb + asymmetry(jg)) … … 379 379 #endif 380 380 381 dojg = 1, ng381 DO jg = 1, ng 382 382 ! Compute upward and downward emission assuming the Planck 383 383 ! function to vary linearly with optical depth within the layer … … 474 474 ! Added for DWD (2020) 475 475 !NEC$ shortloop 476 dojg = 1, ng476 DO jg = 1, ng 477 477 478 478 gamma4 = 1.0_jprd - gamma3(jg) … … 624 624 trans_dir_dir = exp(trans_dir_dir) 625 625 626 dojg = 1, ng626 DO jg = 1, ng 627 627 628 628 ! Zdunkowski "PIFM" (Zdunkowski et al., 1980; Contributions to … … 652 652 exponential = exp(-k_exponent*od) 653 653 654 dojg = 1, ng654 DO jg = 1, ng 655 655 k_mu0 = k_exponent(jg)*mu0 656 656 one_minus_kmu0_sqr = 1.0_jprb - k_mu0*k_mu0 … … 698 698 #else 699 699 ! GPU-capable and vector-optimized version for ICON 700 dojg = 1, ng700 DO jg = 1, ng 701 701 702 702 trans_dir_dir(jg) = max(-max(od(jg) * (1.0_jprb/mu0),0.0_jprb),-1000.0_jprb) … … 810 810 ! Added for DWD (2020) 811 811 !NEC$ shortloop 812 dojg = 1, ng812 DO jg = 1, ng 813 813 ! Note that if the minimum value is reduced (e.g. to 1.0e-24) 814 814 ! then noise starts to appear as a function of solar zenith -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/easy_netcdf.F90
r4773 r5158 392 392 393 393 ndimlens(:) = 0 394 doj = 1,ndims394 DO j = 1,ndims 395 395 istatus = nf90_inquire_dimension(this%ncid, idimids(j), len=ndimlens(j)) 396 396 if (istatus /= NF90_NOERR) then … … 404 404 if (present(ntotal)) then 405 405 ntotal = 1 406 doj = 1, ndims406 DO j = 1, ndims 407 407 ntotal = ntotal * ndimlens(j) 408 408 end do … … 613 613 ! Compute number of elements of the variable in the file 614 614 ntotal = 1 615 doj = 1, ndims615 DO j = 1, ndims 616 616 ntotal = ntotal * ndimlens(j) 617 617 end do … … 658 658 ! Compute number of elements of the variable in the file 659 659 ntotal = 1 660 doj = 1, ndims660 DO j = 1, ndims 661 661 ntotal = ntotal * ndimlens(j) 662 662 end do … … 707 707 ! i.e. excluding the slowest varying dimension, indexed by "index" 708 708 ntotal = 1 709 doj = 1, ndims-1709 DO j = 1, ndims-1 710 710 ntotal = ntotal * ndimlens(j) 711 711 end do … … 761 761 ! Ensure variable has only one dimension in the file 762 762 n = 1 763 doj = 1, ndims763 DO j = 1, ndims 764 764 n = n * ndimlens(j) 765 765 if (j > 1 .and. ndimlens(j) > 1) then … … 819 819 ! Ensure variable has only one dimension in the file 820 820 n = 1 821 doj = 1, ndims821 DO j = 1, ndims 822 822 n = n * ndimlens(j) 823 823 if (j > 1 .and. ndimlens(j) > 1) then … … 878 878 ! Ensure variable has only one dimension in the file 879 879 n = 1 880 doj = 1, ndims880 DO j = 1, ndims 881 881 n = n * ndimlens(j) 882 882 if (j > 1 .and. ndimlens(j) > 1) then … … 938 938 ! Ensure variable has only one dimension aside from the last one 939 939 n = 1 940 doj = 1, ndims-1940 DO j = 1, ndims-1 941 941 n = n * ndimlens(j) 942 942 if (j > 1 .and. ndimlens(j) > 1) then … … 1021 1021 ! dimensions 1022 1022 ntotal = 1 1023 doj = 1, ndims1023 DO j = 1, ndims 1024 1024 ntotal = ntotal * ndimlens(j) 1025 1025 if (j > 2 .and. ndimlens(j) > 1) then … … 1133 1133 ! dimensions 1134 1134 ntotal = 1 1135 doj = 1, ndims1135 DO j = 1, ndims 1136 1136 ntotal = ntotal * ndimlens(j) 1137 1137 if (j > 2 .and. ndimlens(j) > 1) then … … 1202 1202 vcount(1:2) = [ndimlen1,1] 1203 1203 1204 doj = 1,ndimlen21204 DO j = 1,ndimlen2 1205 1205 vstart(2) = j 1206 1206 istatus = nf90_get_var(this%ncid, ivarid, matrix(:,j), start=vstart, count=vcount) … … 1252 1252 ! dimensions aside from the last one 1253 1253 ntotal = 1 1254 doj = 1, ndims-11254 DO j = 1, ndims-1 1255 1255 ntotal = ntotal * ndimlens(j) 1256 1256 if (j > 2 .and. ndimlens(j) > 1) then … … 1376 1376 ! dimensions 1377 1377 ntotal = 1 1378 doj = 1, ndims1378 DO j = 1, ndims 1379 1379 ntotal = ntotal * ndimlens(j) 1380 1380 if (j > 3 .and. ndimlens(j) > 1) then … … 1512 1512 ! dimensions aside from the last one 1513 1513 ntotal = 1 1514 doj = 1, ndims-11514 DO j = 1, ndims-1 1515 1515 ntotal = ntotal * ndimlens(j) 1516 1516 if (j > 3 .and. ndimlens(j) > 1) then … … 1654 1654 ! dimensions 1655 1655 ntotal = 1 1656 doj = 1, ndims1656 DO j = 1, ndims 1657 1657 ntotal = ntotal * ndimlens(j) 1658 1658 if (j > 4 .and. ndimlens(j) > 1) then … … 1804 1804 1805 1805 ! Pad with blanks since nf90_get_att does not do this 1806 doj = i_attr_len+1,len(attr_str)1806 DO j = i_attr_len+1,len(attr_str) 1807 1807 attr_str(j:j) = ' ' 1808 1808 end do … … 1878 1878 1879 1879 ! Pad with blanks since nf90_get_att does not do this 1880 doj = i_attr_len+1,len(attr_str)1880 DO j = i_attr_len+1,len(attr_str) 1881 1881 attr_str(j:j) = ' ' 1882 1882 end do … … 2705 2705 end if 2706 2706 2707 dojdim = 1,ndims2707 DO jdim = 1,ndims 2708 2708 istatus = nf90_inquire_dimension(infile%ncid, idimids(jdim), & 2709 2709 & name=dimname, len=dimlen) … … 2771 2771 2772 2772 ! Map dimension IDs 2773 dojdim = 1,ndims2773 DO jdim = 1,ndims 2774 2774 istatus = nf90_inquire_dimension(infile%ncid, idimids_in(jdim), name=dim_name) 2775 2775 if (istatus /= NF90_NOERR) then … … 2801 2801 2802 2802 ! Copy attributes 2803 dojattr = 1,nattr2803 DO jattr = 1,nattr 2804 2804 istatus = nf90_inq_attname(infile%ncid, ivarid_in, jattr, attr_name) 2805 2805 if (istatus /= NF90_NOERR) then -
LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/utilities/print_matrix.F90
r4773 r5158 16 16 write(unit_local,'(a,a,$)') name, '=[' 17 17 end if 18 doi = 1,size(vec,1)18 DO i = 1,size(vec,1) 19 19 write(unit_local,'(f16.8,$)') vec(i) 20 20 end do … … 43 43 write(unit_local,'(a,a,$)') name, '=[' 44 44 end if 45 doi = 1,size(mat,1)46 doj = 1,size(mat,2)45 DO i = 1,size(mat,1) 46 DO j = 1,size(mat,2) 47 47 write(unit_local,'(f16.8,$)') mat(i,j) 48 48 end do … … 61 61 integer :: i, j, k 62 62 write(6,'(a,a)') name, '=' 63 dok = 1,size(mat,3)64 doi = 1,size(mat,1)65 doj = 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) 66 66 write(6,'(f16.8,$)') mat(i,j,k) 67 67 end do
Note: See TracChangeset
for help on using the changeset viewer.