Changeset 5158 for LMDZ6/branches/Amaury_dev/libf/phylmd/cosp
- Timestamp:
- Aug 2, 2024, 2:12:03 PM (13 months ago)
- Location:
- LMDZ6/branches/Amaury_dev/libf/phylmd/cosp
- Files:
-
- 30 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/array_lib.F90
r5099 r5158 141 141 xsort = xarr(ist) 142 142 143 doi=1,nxx143 DO i=1,nxx 144 144 iloc = infind(xsort,xxarr(i),dist=d) 145 145 if (d > tol) then -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/congvec.h
r5099 r5158 37 37 ! *****************************COPYRIGHT******************************* 38 38 39 doirand = 1, npoints39 DO irand = 1, npoints 40 40 ! Marsaglia CONG algorithm 41 41 seed(irand)=69069*seed(irand)+1234567 … … 48 48 overflow_32=i2_16*i2_16 49 49 if ( overflow_32 .le. huge32 ) then 50 doirand = 1, npoints50 DO irand = 1, npoints 51 51 ran(irand)=ran(irand)+1 52 52 ran(irand)=(ran(irand))-int(ran(irand)) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/cosp_output_mod.F90
r5157 r5158 72 72 "clmcalipsoice", "CALIPSO Ice-Phase Mid Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 73 73 TYPE(ctrl_outcosp), SAVE :: o_clmcalipsoliq = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 74 "clmcalipsoliq", "CALIPSO Liq-Phase Mid Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 74 "clmcalipsoliq", "CALIPSO Liq-Phase Mid Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 75 75 TYPE(ctrl_outcosp), SAVE :: o_clhcalipsoice = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 76 76 "clhcalipsoice", "CALIPSO Ice-Phase High Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) … … 80 80 "cltcalipsoice", "CALIPSO Ice-Phase Tot Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 81 81 TYPE(ctrl_outcosp), SAVE :: o_cltcalipsoliq = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 82 "cltcalipsoliq", "CALIPSO Liq-Phase Tot Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 82 "cltcalipsoliq", "CALIPSO Liq-Phase Tot Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) 83 83 TYPE(ctrl_outcosp), SAVE :: o_cllcalipsoun = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 84 84 "cllcalipsoun", "CALIPSO Undefined-Phase Low Level Cloud Fraction", "%", (/ ('', i=1, 3) /)) … … 94 94 "clcalipsoliq", "Lidar Liq-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /)) 95 95 TYPE(ctrl_outcosp), SAVE :: o_clcalipsoun = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 96 "clcalipsoun", "Lidar Undef-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /)) 96 "clcalipsoun", "Lidar Undef-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /)) 97 97 TYPE(ctrl_outcosp), SAVE :: o_clcalipsotmpice = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & 98 98 "clcalipsotmpice", "Lidar Ice-Phase Cloud Fraction", "%", (/ ('', i=1, 3) /)) … … 115 115 "clcalipsothin", "Lidar Thin profile Cloud Fraction", "%", (/ ('', i=1, 3) /)) !OPAQ 116 116 TYPE(ctrl_outcosp), SAVE :: o_clcalipsozopaque = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & !OPAQ 117 "clcalipsozopaque", "Lidar z_opaque Fraction", "%", (/ ('', i=1, 3) /)) 117 "clcalipsozopaque", "Lidar z_opaque Fraction", "%", (/ ('', i=1, 3) /)) !OPAQ 118 118 TYPE(ctrl_outcosp), SAVE :: o_clcalipsoopacity = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & !OPAQ 119 "clcalipsoopacity", "Lidar opacity Fraction", "%", (/ ('', i=1, 3) /)) 119 "clcalipsoopacity", "Lidar opacity Fraction", "%", (/ ('', i=1, 3) /)) !OPAQ 120 120 121 121 TYPE(ctrl_outcosp), SAVE :: o_proftemp = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), & !TIBO … … 243 243 integer :: Nlevlmdz, Ncolumns ! Number of levels 244 244 real,dimension(Nlevlmdz) :: presnivs 245 real :: dtime, freq_cosp, ecrit_day, ecrit_hf, ecrit_mth245 REAL :: dtime, freq_cosp, ecrit_day, ecrit_hf, ecrit_mth 246 246 logical :: ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, use_vgrid, ok_all_xml 247 247 type(cosp_vgrid) :: vgrid ! Information on vertical grid of stats … … 250 250 !!! Variables locales 251 251 integer :: idayref, iff, ii 252 real:: zjulian,zjulian_start252 REAL :: zjulian,zjulian_start 253 253 real,dimension(Ncolumns) :: column_ax 254 254 real,dimension(DBZE_BINS) :: dbze_ax … … 271 271 272 272 !! Definition valeurs axes 273 doii=1,Ncolumns273 DO ii=1,Ncolumns 274 274 column_ax(ii) = real(ii) 275 275 enddo 276 276 277 doi=1,DBZE_BINS277 DO i=1,DBZE_BINS 278 278 dbze_ax(i) = CFAD_ZE_MIN + CFAD_ZE_WIDTH*(i - 0.5) 279 279 enddo -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/cosp_output_write_mod.F90
r5157 r5158 31 31 !!! Variables d'entree 32 32 integer :: itap, Nlevlmdz, Ncolumns, Npoints 33 real:: freq_COSP, dtime, missing_val, missing_cosp33 REAL :: freq_COSP, dtime, missing_val, missing_cosp 34 34 type(cosp_config) :: cfg ! Control outputs 35 35 type(cosp_gridbox) :: gbx ! Gridbox information. Input for COSP … … 161 161 162 162 IF (using_xios) THEN 163 doicl=1,SR_BINS163 DO icl=1,SR_BINS 164 164 tmp_fi4da_cfadL(:,:,icl)=stlidar%cfad_sr(:,icl,:) 165 165 enddo … … 169 169 ELSE 170 170 if (cfg%LcfadLidarsr532) then 171 doicl=1,SR_BINS171 DO icl=1,SR_BINS 172 172 CALL histwrite3d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr(:,icl,:),nvert,icl) 173 173 enddo 174 174 endif 175 175 if (cfg%LprofSR) then 176 doicl=1,Ncolumns !TIBO176 DO icl=1,Ncolumns !TIBO 177 177 CALL histwrite3d_cosp(o_profSR,stlidar%profSR(:,icl,:),nvert,icl) !TIBO 178 178 enddo !TIBO … … 183 183 184 184 if (cfg%LparasolRefl) then 185 dok=1,PARASOL_NREFL186 doip=1, Npoints185 DO k=1,PARASOL_NREFL 186 DO ip=1, Npoints 187 187 if (stlidar%cldlayer(ip,4).gt.0.01.and.stlidar%parasolrefl(ip,k).ne.missing_val) then 188 188 parasolcrefl(ip,k)=(stlidar%parasolrefl(ip,k)-0.03*(1.-stlidar%cldlayer(ip,4)))/ & … … 203 203 ELSE 204 204 if (cfg%Latb532) then 205 doicl=1,Ncolumns205 DO icl=1,Ncolumns 206 206 CALL histwrite3d_cosp(o_atb532,sglidar%beta_tot(:,icl,:),nvertmcosp,icl) 207 207 enddo … … 218 218 where(stradar%cfad_ze == R_UNDEF) stradar%cfad_ze = missing_val 219 219 IF (using_xios) THEN 220 doicl=1,DBZE_BINS220 DO icl=1,DBZE_BINS 221 221 tmp_fi4da_cfadR(:,:,icl)=stradar%cfad_ze(:,icl,:) 222 222 enddo … … 226 226 ELSE 227 227 if (cfg%Ldbze94) then 228 doicl=1,Ncolumns228 DO icl=1,Ncolumns 229 229 CALL histwrite3d_cosp(o_dbze94,sgradar%Ze_tot(:,icl,:),nvert,icl) 230 230 enddo 231 231 endif 232 232 if (cfg%LcfadDbze94) then 233 doicl=1,DBZE_BINS233 DO icl=1,DBZE_BINS 234 234 CALL histwrite3d_cosp(o_cfadDbze94,stradar%cfad_ze(:,icl,:),nvert,icl) 235 235 enddo … … 266 266 ELSE 267 267 if (cfg%Lclisccp) then 268 doicl=1,7268 DO icl=1,7 269 269 CALL histwrite3d_cosp(o_clisccp2,isccp%fq_isccp(:,icl,:),nvertisccp,icl) 270 270 enddo … … 287 287 288 288 IF (using_xios) THEN 289 doicl=1,MISR_N_CTH289 DO icl=1,MISR_N_CTH 290 290 tmp_fi4da_misr(:,icl,:)=misr%fq_MISR(:,:,icl) 291 291 enddo … … 294 294 ELSE 295 295 if (cfg%LclMISR) then 296 doicl=1,7296 DO icl=1,7 297 297 CALL histwrite3d_cosp(o_clMISR,misr%fq_MISR(:,icl,:),nvertmisr,icl) 298 298 enddo … … 367 367 ELSE 368 368 if (cfg%Lclmodis) then 369 doicl=1,7369 DO icl=1,7 370 370 CALL histwrite3d_cosp(o_clmodis, & 371 371 modis%Optical_Thickness_vs_Cloud_Top_Pressure(:,icl,:),nvertisccp,icl) … … 385 385 ELSE 386 386 if (cfg%Lcrimodis) then 387 doicl=1,7387 DO icl=1,7 388 388 CALL histwrite3d_cosp(o_crimodis, & 389 389 modis%Optical_Thickness_vs_ReffIce(:,icl,:),nvertReffIce,icl) … … 391 391 endif 392 392 if (cfg%Lcrlmodis) then 393 doicl=1,7393 DO icl=1,7 394 394 CALL histwrite3d_cosp(o_crlmodis, & 395 395 modis%Optical_Thickness_vs_ReffLiq(:,icl,:),nvertReffLiq,icl) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/cosp_read_otputkeys.F90
r5133 r5158 28 28 29 29 30 doi=1,N_OUT_LIST30 DO i=1,N_OUT_LIST 31 31 cfg%out_list(i)='' 32 32 enddo … … 130 130 131 131 132 doi=1,N_OUT_LIST132 DO i=1,N_OUT_LIST 133 133 cfg%out_list(i)='' 134 134 enddo … … 269 269 LprofSR,Lproftemp !TIBO (2) 270 270 271 doi=1,N_OUT_LIST271 DO i=1,N_OUT_LIST 272 272 cfg%out_list(i)='' 273 273 enddo … … 771 771 IF (using_xios) THEN 772 772 773 doi=1,N_OUT_LIST773 DO i=1,N_OUT_LIST 774 774 cfg%out_list(i)='' 775 775 enddo -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/dsd.F90
r5099 r5158 311 311 lidx = infind(D,dmin) 312 312 uidx = infind(D,dmax) 313 dok=lidx,uidx313 DO k=lidx,uidx 314 314 315 315 N(k) = ( & -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/format_input.F90
r5099 r5158 52 52 53 53 nprof = size(hgt_matrix,1) 54 doi=1,nprof54 DO i=1,nprof 55 55 call lin_interpolate(env_t_matrix(i,:),env_hgt_matrix(i,:), & 56 56 t_matrix(i,:),hgt_matrix(i,:),1000._KR8) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/gases.F90
r5099 r5158 126 126 sumo = 0. 127 127 aux1 = 1.1*e_th 128 doi=1,nbands_o2128 DO i=1,nbands_o2 129 129 aux2 = f/v0(i) 130 130 aux3 = v0(i)-f … … 151 151 sumo = 0. 152 152 aux1 = 4.8*e_th 153 doi=1,nbands_h2o153 DO i=1,nbands_h2o 154 154 aux2 = f/v1(i) 155 155 aux3 = v1(i)-f -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/lidar_simulator.F90
r5099 r5158 324 324 ! altitude at half pressure levels: 325 325 zheight(:,1) = 0.0 326 dok = 2, nlev+1326 DO k = 2, nlev+1 327 327 zheight(:,k) = zheight(:,k-1) & 328 328 -(presf(:,k)-presf(:,k-1))/(rhoair(:,k-1)*9.81) … … 341 341 342 342 ! polynomes kp_lidar derived from Mie theory: 343 doi = 1, npart343 DO i = 1, npart 344 344 where ( rad_part(:,:,i).gt.0.0) 345 345 kp_part(:,:,i) = & … … 361 361 362 362 ! alpha of particles in each subcolumn: 363 doi = 1, npart363 DO i = 1, npart 364 364 where ( rad_part(:,:,i).gt.0.0) 365 365 alpha_part(:,:,i) = 3.0/4.0 * Qscat & -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/math_lib.F90
r5099 r5158 49 49 ga=1.0 50 50 m1=x-1 51 dok=2,m151 DO k=2,m1 52 52 ga=ga*k 53 53 enddo … … 60 60 m=int(z) 61 61 r=1.0 62 dok=1,m62 DO k=1,m 63 63 r=r*(z-k) 64 64 enddo … … 80 80 -.206d-13, -.54d-14, .14d-14, .1d-15/ 81 81 gr=g(26) 82 dok=25,1,-182 DO k=25,1,-1 83 83 gr=gr*z+g(k) 84 84 enddo … … 146 146 else 147 147 sumo = 0. 148 doj=i1,i2148 DO j=i1,i2 149 149 deltah = abs(s(i1+1)-s(i1)) 150 150 sumo = sumo + f(j)*deltah … … 268 268 end if 269 269 270 doi = 2, ntab270 DO i = 2, ntab 271 271 if ( xtab(i) <= xtab(i-1) ) then 272 272 lerror = .true. … … 311 311 ihi = ntab 312 312 313 doi = 1, ntab313 DO i = 1, ntab 314 314 if ( a <= xtab(i) ) then 315 315 exit … … 321 321 ilo = min ( ilo, ntab - 1 ) 322 322 323 doi = 1, ntab323 DO i = 1, ntab 324 324 if ( xtab(i) <= b ) then 325 325 exit … … 335 335 sum1 = 0.0D+00 336 336 337 doi = ilo, ihi337 DO i = ilo, ihi 338 338 339 339 x1 = xtab(i-1) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp.F90
r5099 r5158 68 68 integer,dimension(2) :: ix,iy 69 69 logical :: reff_zero 70 real:: maxp,minp70 REAL :: maxp,minp 71 71 integer,dimension(:),allocatable :: & ! Dimensions nPoints 72 72 seed ! It is recommended that the seed is set to a different value for each model … … 95 95 96 96 !++++++++++ Depth of model layers ++++++++++++ 97 doi=1,Nlevels-197 DO i=1,Nlevels-1 98 98 gbx%dlev(:,i) = gbx%zlev_half(:,i+1) - gbx%zlev_half(:,i) 99 99 enddo … … 195 195 Niter = gbx%Npoints/gbx%Npoints_it ! Integer division 196 196 if (Niter*gbx%Npoints_it < gbx%Npoints) Niter = Niter + 1 197 doi=1,Niter197 DO i=1,Niter 198 198 i_first = (i-1)*gbx%Npoints_it + 1 199 199 i_last = i_first + gbx%Npoints_it - 1 … … 409 409 410 410 ! Precipitation fraction 411 doj=1,Npoints,1412 dok=1,Nlevels,1413 doi=1,Ncolumns,1411 DO j=1,Npoints,1 412 DO k=1,Nlevels,1 413 DO i=1,Ncolumns,1 414 414 if (sgx%frac_out (j,i,Nlevels+1-k) == I_LSC) frac_ls(j,k)=frac_ls(j,k)+1. 415 415 if (sgx%frac_out (j,i,Nlevels+1-k) == I_CVC) frac_cv(j,k)=frac_cv(j,k)+1. … … 434 434 else 435 435 ! This is done within a loop (unvectorized) over nPoints to save memory 436 doj=1,Npoints436 DO j=1,Npoints 437 437 sgx%frac_out(j,:,1:Nlevels) = sgx%frac_out(j,:,Nlevels:1:-1) 438 438 sgx%prec_frac(j,:,1:Nlevels) = sgx%prec_frac(j,:,Nlevels:1:-1) … … 445 445 ! Populate the subgrid arrays 446 446 call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,Nhydro,sghydro) 447 dok=1,Ncolumns447 DO k=1,Ncolumns 448 448 !--------- Mixing ratios for clouds and Reff for Clouds and precip ------- 449 449 column_frac_out => sgx%frac_out(:,k,:) … … 498 498 enddo 499 499 ! convert the mixing ratio and precipitation flux from gridbox mean to the fraction-based values 500 dok=1,Nlevels501 doj=1,Npoints500 DO k=1,Nlevels 501 DO j=1,Npoints 502 502 !--------- Clouds ------- 503 503 if (frac_ls(j,k) .ne. 0.) then -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_constants.F90
r5099 r5158 164 164 integer :: HCLASS_TYPE(N_HYDRO),HCLASS_PHASE(N_HYDRO) 165 165 166 real:: HCLASS_DMIN(N_HYDRO),HCLASS_DMAX(N_HYDRO), &166 REAL :: HCLASS_DMIN(N_HYDRO),HCLASS_DMAX(N_HYDRO), & 167 167 HCLASS_APM(N_HYDRO),HCLASS_BPM(N_HYDRO),HCLASS_RHO(N_HYDRO), & 168 168 HCLASS_P1(N_HYDRO),HCLASS_P2(N_HYDRO),HCLASS_P3(N_HYDRO) … … 274 274 integer :: HCLASS_TYPE(N_HYDRO),HCLASS_PHASE(N_HYDRO) 275 275 276 real :: HCLASS_DMIN(N_HYDRO),HCLASS_DMAX(N_HYDRO), &276 REAL :: HCLASS_DMIN(N_HYDRO),HCLASS_DMAX(N_HYDRO), & 277 277 HCLASS_APM(N_HYDRO),HCLASS_BPM(N_HYDRO),HCLASS_RHO(N_HYDRO), & 278 278 HCLASS_P1(N_HYDRO),HCLASS_P2(N_HYDRO),HCLASS_P3(N_HYDRO) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_isccp_simulator.F90
r5099 r5158 45 45 ! Local variables 46 46 integer :: Nlevels,Npoints 47 real:: pfull(gbx%Npoints, gbx%Nlevels)48 real:: phalf(gbx%Npoints, gbx%Nlevels + 1)49 real:: qv(gbx%Npoints, gbx%Nlevels)50 real:: cc(gbx%Npoints, gbx%Nlevels)51 real:: conv(gbx%Npoints, gbx%Nlevels)52 real:: dtau_s(gbx%Npoints, gbx%Nlevels)53 real:: dtau_c(gbx%Npoints, gbx%Nlevels)54 real:: at(gbx%Npoints, gbx%Nlevels)55 real:: dem_s(gbx%Npoints, gbx%Nlevels)56 real:: dem_c(gbx%Npoints, gbx%Nlevels)57 real:: frac_out(gbx%Npoints, gbx%Ncolumns, gbx%Nlevels)47 REAL :: pfull(gbx%Npoints, gbx%Nlevels) 48 REAL :: phalf(gbx%Npoints, gbx%Nlevels + 1) 49 REAL :: qv(gbx%Npoints, gbx%Nlevels) 50 REAL :: cc(gbx%Npoints, gbx%Nlevels) 51 REAL :: conv(gbx%Npoints, gbx%Nlevels) 52 REAL :: dtau_s(gbx%Npoints, gbx%Nlevels) 53 REAL :: dtau_c(gbx%Npoints, gbx%Nlevels) 54 REAL :: at(gbx%Npoints, gbx%Nlevels) 55 REAL :: dem_s(gbx%Npoints, gbx%Nlevels) 56 REAL :: dem_c(gbx%Npoints, gbx%Nlevels) 57 REAL :: frac_out(gbx%Npoints, gbx%Ncolumns, gbx%Nlevels) 58 58 integer :: sunlit(gbx%Npoints) 59 59 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_lidar.F90
r5099 r5158 54 54 ! Local variables 55 55 integer :: i 56 real:: presf(sgx%Npoints, sgx%Nlevels + 1)56 REAL :: presf(sgx%Npoints, sgx%Nlevels + 1) 57 57 real,dimension(sgx%Npoints, sgx%Nlevels) :: lsca,mr_ll,mr_li,mr_cl,mr_ci 58 58 real,dimension(sgx%Npoints, sgx%Nlevels) :: beta_tot,tau_tot … … 63 63 presf(:,sgx%Nlevels + 1) = 0.0 64 64 lsca = gbx%tca-gbx%cca 65 doi=1,sgx%Ncolumns65 DO i=1,sgx%Ncolumns 66 66 ! Temporary arrays for simulator call 67 67 mr_ll(:,:) = sghydro%mr_hydro(:,i,:,I_LSCLIQ) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_misr_simulator.F90
r5099 r5158 49 49 ! Local variables 50 50 integer :: Nlevels,Npoints 51 real:: dtau_s(gbx%Npoints, gbx%Nlevels)52 real:: dtau_c(gbx%Npoints, gbx%Nlevels)53 real:: at(gbx%Npoints, gbx%Nlevels)54 real:: frac_out(gbx%Npoints, gbx%Ncolumns, gbx%Nlevels)51 REAL :: dtau_s(gbx%Npoints, gbx%Nlevels) 52 REAL :: dtau_c(gbx%Npoints, gbx%Nlevels) 53 REAL :: at(gbx%Npoints, gbx%Nlevels) 54 REAL :: frac_out(gbx%Npoints, gbx%Ncolumns, gbx%Nlevels) 55 55 integer :: sunlit(gbx%Npoints) 56 56 57 real:: zfull(gbx%Npoints, gbx%Nlevels) ! height (in meters) of full model levels (i.e. midpoints)57 REAL :: zfull(gbx%Npoints, gbx%Nlevels) ! height (in meters) of full model levels (i.e. midpoints) 58 58 ! zfull(npoints,1) is top level of model 59 59 ! zfull(npoints,nlev) is bottom level of model -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_modis_simulator.F90
r5099 r5158 162 162 163 163 ! Loop version of spread above - intrinsic doesn't work on certain platforms. 164 dok = 1, nLevels165 doj = 1, nSubCols166 doi = 1, nSunlit164 DO k = 1, nLevels 165 DO j = 1, nSubCols 166 DO i = 1, nSunlit 167 167 if(subCols%frac_out(sunlit(i), j, k) == I_LSC) then 168 168 opticalThickness(i, j, k) = gridBox%dtau_s(sunlit(i), k) … … 186 186 187 187 ! Loop version of spread above - intrinsic doesn't work on certain platforms. 188 dok = 1, nLevels189 doj = 1, nSubCols190 doi = 1, nSunlit188 DO k = 1, nLevels 189 DO j = 1, nSubCols 190 DO i = 1, nSunlit 191 191 if(subCols%frac_out(sunlit(i), j, k) == I_CVC) opticalThickness(i, j, k) = gridBox%dtau_c(sunlit(i), k) 192 192 end do … … 205 205 isccpCloudTopPressure(:, :) = isccpSim%boxptop(sunlit(:), :) 206 206 207 doi = 1, nSunlit207 DO i = 1, nSunlit 208 208 call modis_L2_simulator(temperature(i, :), pressureLayers(i, :), pressureLevels(i, :), & 209 209 opticalThickness(i, :, :), cloudWater(i, :, :), cloudIce(i, :, :), & -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_radar.F90
r5099 r5158 139 139 140 140 ! set flag denoting position of radar relative to hgt_matrix orientation 141 141 ngate = size(hgt_matrix,2) 142 142 143 143 hgt_descending = hgt_matrix(1,1) > hgt_matrix(1,ngate) 144 144 145 146 147 148 149 150 151 152 153 145 if ( & 146 (gbx%surface_radar == 1 .and. hgt_descending) .or. & 147 (gbx%surface_radar == 0 .and. (.not. hgt_descending)) & 148 ) & 149 then 150 gbx%hp%radar_at_layer_one = .false. 151 else 152 gbx%hp%radar_at_layer_one = .true. 153 endif 154 154 155 155 ! ----- loop over subcolumns ----- 156 dopr=1,sgx%Ncolumns156 DO pr=1,sgx%Ncolumns 157 157 158 158 ! NOTE: … … 160 160 ! only hydrometeor profiles will be different for each subgridbox 161 161 162 doi=1,gbx%Nhydro162 DO i=1,gbx%Nhydro 163 163 hm_matrix(i,:,:) = sghydro%mr_hydro(:,pr,:,i)*1000.0 ! Units from kg/kg to g/kg 164 164 if (gbx%use_reff) then -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_stats.F90
r5099 r5158 203 203 logical :: lunits 204 204 integer :: l 205 real:: w ! Weight206 real:: dbb, dtb, dbt, dtt ! Distances between edges of both grids205 REAL :: w ! Weight 206 REAL :: dbb, dtb, dbt, dtt ! Distances between edges of both grids 207 207 integer :: Nw ! Number of weights 208 real:: wt ! Sum of weights208 REAL :: wt ! Sum of weights 209 209 real,dimension(Nlevels) :: oldgrid_bot,oldgrid_top ! Lower and upper boundaries of model grid 210 real:: yp ! Local copy of y at a particular point.210 REAL :: yp ! Local copy of y at a particular point. 211 211 ! This allows for change of units. 212 212 … … 216 216 r = 0.0 217 217 218 doi=1,Npoints218 DO i=1,Npoints 219 219 ! Calculate tops and bottoms of new and old grids 220 220 oldgrid_bot = zhalf(i,:) … … 223 223 l = 0 ! Index of level in the old grid 224 224 ! Loop over levels in the new grid 225 dok = 1,Nglevels225 DO k = 1,Nglevels 226 226 Nw = 0 ! Number of weigths 227 227 wt = 0.0 ! Sum of weights 228 228 ! Loop over levels in the old grid and accumulate total for weighted average 229 do229 DO 230 230 l = l + 1 231 231 w = 0.0 ! Initialise weight to 0 … … 254 254 Nw = Nw + 1 255 255 wt = wt + w 256 doj=1,Ncolumns256 DO j=1,Ncolumns 257 257 if (lunits) then 258 258 if (y(i,j,l) /= R_UNDEF) then … … 273 273 ! Calculate average in new grid 274 274 if (Nw > 0) then 275 doj=1,Ncolumns275 DO j=1,Ncolumns 276 276 r(i,j,k) = r(i,j,k)/wt 277 277 enddo … … 281 281 282 282 ! Set points under surface to R_UNDEF, and change to dBZ if necessary 283 dok=1,Nglevels284 doj=1,Ncolumns285 doi=1,Npoints283 DO k=1,Nglevels 284 DO j=1,Ncolumns 285 DO i=1,Npoints 286 286 if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level 287 287 if (lunits) then -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_types.F90
r5099 r5158 44 44 Lclcalipsoliq,Lclcalipsoice,Lclcalipsoun, & 45 45 Lclcalipsotmp,Lclcalipsotmpliq,Lclcalipsotmpice,Lclcalipsotmpun, & 46 46 Lcltcalipsoliq,Lcltcalipsoice,Lcltcalipsoun, & 47 47 Lclhcalipsoliq,Lclhcalipsoice,Lclhcalipsoun, & 48 48 Lclmcalipsoliq,Lclmcalipsoice,Lclmcalipsoun, & … … 274 274 275 275 ! Radar ancillary info 276 real:: radar_freq, & ! Radar frequency [GHz]276 REAL :: radar_freq, & ! Radar frequency [GHz] 277 277 k2 ! |K|^2, -1=use frequency dependent default 278 278 integer :: surface_radar, & ! surface=1, spaceborne=0 … … 383 383 ! 2 = experimental setting 384 384 integer :: isccp_overlap ! overlap type (1=max, 2=rand, 3=max/rand) 385 real:: isccp_emsfc_lw ! 10.5 micron emissivity of surface (fraction)385 REAL :: isccp_emsfc_lw ! 10.5 micron emissivity of surface (fraction) 386 386 387 387 ! RTTOV inputs/options … … 392 392 integer, dimension(:), pointer :: Ichan ! Channel numbers 393 393 real, dimension(:), pointer :: Surfem ! Surface emissivity 394 real:: ZenAng ! Satellite Zenith Angles395 real:: co2,ch4,n2o,co ! Mixing ratios of trace gases394 REAL :: ZenAng ! Satellite Zenith Angles 395 REAL :: co2,ch4,n2o,co ! Mixing ratios of trace gases 396 396 397 397 END TYPE COSP_GRIDBOX … … 550 550 ! Local variables 551 551 integer :: i 552 real:: zstep552 REAL :: zstep 553 553 554 554 x%use_vgrid = use_vgrid … … 587 587 zstep = 20000.0/x%Nlvgrid 588 588 endif 589 doi=1,x%Nlvgrid589 DO i=1,x%Nlvgrid 590 590 x%zl(i) = (i-1)*zstep 591 591 x%zu(i) = i*zstep … … 654 654 x%tau_tot = 0.0 655 655 x%refl = 0.0 ! parasol 656 x%temp_tot 657 x%betaperp_tot = 0.0656 x%temp_tot = 0.0 657 x%betaperp_tot = 0.0 658 658 END SUBROUTINE CONSTRUCT_COSP_SGLIDAR 659 659 … … 1204 1204 ! y%hp%idd = x%hp%idd 1205 1205 sz = shape(x%hp%Z_scale_flag) 1206 dok=1,sz(3)1207 doj=1,sz(2)1208 doi=1,sz(1)1206 DO k=1,sz(3) 1207 DO j=1,sz(2) 1208 DO i=1,sz(1) 1209 1209 if (x%hp%N_scale_flag(i,k)) y%hp%N_scale_flag(i,k) = .true. 1210 1210 if (x%hp%Z_scale_flag(i,j,k)) y%hp%Z_scale_flag(i,j,k) = .true. -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_utils.F90
r5099 r5158 58 58 ! Local variables 59 59 integer :: i,j,k 60 real:: sigma,one_over_xip1,xi,rho0,rho,lambda_x,gamma_4_3_2,delta61 real :: seuil60 REAL :: sigma,one_over_xip1,xi,rho0,rho,lambda_x,gamma_4_3_2,delta 61 REAL :: seuil 62 62 63 63 if (ok_debug_cosp) then … … 77 77 delta = (alpha_x + b_x + d_x - n_bx + 1.0) 78 78 79 dok=1,Nlevels80 doj=1,Ncolumns81 doi=1,Npoints79 DO k=1,Nlevels 80 DO j=1,Ncolumns 81 DO i=1,Npoints 82 82 if ((prec_frac(i,j,k)==prec_type).or.(prec_frac(i,j,k)==3.)) then 83 83 rho = p(i,k)/(287.05*T(i,k)) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_llnl_stats.F90
r5099 r5158 67 67 ! Valid x values smaller than bmin and larger than bmax are set 68 68 ! into the smallest bin and largest bin, respectively. 69 doj = 1, Nlevels, 170 dok = 1, Ncolumns, 171 doi = 1, Npoints, 169 DO j = 1, Nlevels, 1 70 DO k = 1, Ncolumns, 1 71 DO i = 1, Npoints, 1 72 72 if (x(i,k,j) == R_GROUND) then 73 73 cosp_cfad(i,:,j) = R_UNDEF … … 101 101 102 102 ! local variables 103 real:: sc_ratio104 real:: s_cld, s_att103 REAL :: sc_ratio 104 REAL :: s_cld, s_att 105 105 parameter (S_cld = 5.0) 106 106 parameter (s_att = 0.01) … … 111 111 lidar_only_freq_cloud = 0.0 112 112 tcc = 0.0 113 dopr=1,Npoints114 doi=1,Ncolumns113 DO pr=1,Npoints 114 DO i=1,Ncolumns 115 115 flag_sat = 0 116 116 flag_cld = 0 117 doj=Nlevels,1,-1 !top->surf117 DO j=Nlevels,1,-1 !top->surf 118 118 sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j) 119 119 if ((sc_ratio .le. s_att) .and. (flag_sat .eq. 0)) flag_sat = j -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_lmd_ipsl_stats.F90
r5099 r5158 131 131 132 132 ! SR detection threshold 133 real, parameter :: S_cld_att = 30. ! New threshold for undefine cloud phase detection 133 real, parameter :: S_cld_att = 30. ! New threshold for undefine cloud phase detection 134 134 135 135 ! c ------------------------------------------------------- … … 144 144 ! c ------------------------------------------------------- 145 145 146 doic = 1, ncol146 DO ic = 1, ncol 147 147 pnorm_c = pnorm(:,ic,:) 148 148 where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) … … 152 152 end where 153 153 x3d(:,ic,:) = x3d_c 154 ! 155 154 ! profSR(:,ic,:) = x3d(:,ic,:) !TIBO 155 profSR(:,:,ic) = x3d(:,ic,:) !TIBO2 156 156 enddo 157 157 … … 188 188 parasolrefl(:,:) = 0.0 189 189 190 dok = 1, nrefl191 doic = 1, ncol190 DO k = 1, nrefl 191 DO ic = 1, ncol 192 192 parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k) 193 193 enddo 194 194 enddo 195 195 196 dok = 1, nrefl196 DO k = 1, nrefl 197 197 parasolrefl(:,k) = parasolrefl(:,k) / float(ncol) 198 198 ! if land=1 -> parasolrefl=undef … … 253 253 srbval(5) = 7.0 254 254 srbval(6) = 10.0 255 doi = 7, MIN(10,Nbins)255 DO i = 7, MIN(10,Nbins) 256 256 srbval(i) = srbval(i-1) + 5.0 257 257 enddo … … 268 268 ! c c- Compute CFAD 269 269 ! c ------------------------------------------------------- 270 doj = 1, Nlevels271 doib = 1, Nbins272 dok = 1, Ncolumns273 doi = 1, Npoints270 DO j = 1, Nlevels 271 DO ib = 1, Nbins 272 DO k = 1, Ncolumns 273 DO i = 1, Npoints 274 274 if (x(i,k,j) /= undef) then 275 275 if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) & … … 310 310 integer,parameter :: Ntemp=40 ! indice of the temperature vector 311 311 integer ip, k, iz, ic, ncol, nlev, i, itemp ! loop indice 312 real S_cld_att ! New threshold for undefine cloud phase detection (SR=30) 312 real S_cld_att ! New threshold for undefine cloud phase detection (SR=30) 313 313 integer toplvlsat ! level of the first cloud with SR>30 314 314 real alpha50, beta50, gamma50, delta50, epsilon50, zeta50 ! Polynomial Coef of the phase … … 316 316 317 317 ! Input variables 318 real tmp(Npoints,Nlevels) 318 real tmp(Npoints,Nlevels) ! temperature 319 319 real ATB(Npoints,Ncolumns,Nlevels) ! 3D Attenuated backscatter 320 320 real ATBperp(Npoints,Ncolumns,Nlevels) ! 3D perpendicular attenuated backscatter … … 332 332 333 333 ! Local variables 334 real tmpi(Npoints,Ncolumns,Nlevels) 335 real tmpl(Npoints,Ncolumns,Nlevels) 336 real tmpu(Npoints,Ncolumns,Nlevels) 334 real tmpi(Npoints,Ncolumns,Nlevels) ! temperature of ice cld 335 real tmpl(Npoints,Ncolumns,Nlevels) ! temperature of liquid cld 336 real tmpu(Npoints,Ncolumns,Nlevels) ! temperature of undef cld 337 337 338 338 real checktemp, ATBperp_tmp ! temporary variable … … 399 399 -54.,-51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18., & 400 400 -15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,21.,24.,200. /) 401 401 402 402 ! convert C to K 403 403 tempmod=tempmod+273.15 … … 418 418 ! --------------------------------------------------------------- 419 419 420 dok = 1, Nlevels420 DO k = 1, Nlevels 421 421 422 422 ! cloud detection at subgrid-scale: … … 454 454 nsublay(:,:,4) = 0.0 455 455 456 dok = Nlevels, 1, -1457 doic = 1, Ncolumns458 doip = 1, Npoints456 DO k = Nlevels, 1, -1 457 DO ic = 1, Ncolumns 458 DO ip = 1, Npoints 459 459 460 460 if(srok(ip,ic,k).gt.0.)then 461 461 ! Computation of the cloud fraction as a function of the temperature 462 462 ! instead of height, for ice,liquid and all clouds 463 doitemp=1,Ntemp463 DO itemp=1,Ntemp 464 464 if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then 465 465 lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1. … … 469 469 470 470 if (cldy(ip,ic,k).eq.1.) then 471 doitemp=1,Ntemp471 DO itemp=1,Ntemp 472 472 if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then 473 473 lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1. … … 505 505 cldlay = 0.0 506 506 nsublay = 0.0 507 dok = Nlevels, 1, -1508 doic = 1, Ncolumns509 doip = 1, Npoints507 DO k = Nlevels, 1, -1 508 DO ic = 1, Ncolumns 509 DO ip = 1, Npoints 510 510 511 511 ! Computation of the cloud fraction as a function of the temperature 512 512 ! instead of height, for ice,liquid and all clouds 513 513 if(srok(ip,ic,k).gt.0.)then 514 doitemp=1,Ntemp514 DO itemp=1,Ntemp 515 515 if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then 516 516 lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1. … … 520 520 521 521 if(cldy(ip,ic,k).eq.1.)then 522 doitemp=1,Ntemp522 DO itemp=1,Ntemp 523 523 if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then 524 524 lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1. … … 562 562 nsublayer = 0.0 563 563 564 doiz = 1, Ncat565 doic = 1, Ncolumns564 DO iz = 1, Ncat 565 DO ic = 1, Ncolumns 566 566 567 567 cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz) … … 582 582 ! 4.1 - For Cloudy pixels with 8.16km < z < 19.2km 583 583 ! --------------------------------------------------------------- 584 doncol=1,Ncolumns585 doi=1,Npoints586 587 donlev=Nlevels,18,-1 ! from 19.2km until 8.16km584 DO ncol=1,Ncolumns 585 DO i=1,Npoints 586 587 DO nlev=Nlevels,18,-1 ! from 19.2km until 8.16km 588 588 p1 = pplay(i,nlev) 589 589 590 590 591 591 ! Avoid zero values 592 592 if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then 593 593 ! Computation of the ATBperp along the phase discrimination line 594 594 ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + & … … 609 609 lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1. ! keep the information "temperature criterium used" 610 610 ! to classify the phase cloud 611 611 cldlayphase(i,ncol,4,2) = 1. ! tot cloud 612 612 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 613 614 615 616 617 613 cldlayphase(i,ncol,3,2) = 1. 614 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 615 cldlayphase(i,ncol,2,2) = 1. 616 else ! low cloud 617 cldlayphase(i,ncol,1,2) = 1. 618 618 endif 619 620 621 622 623 624 625 619 cldlayphase(i,ncol,4,5) = 1. ! tot cloud 620 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 621 cldlayphase(i,ncol,3,5) = 1. 622 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 623 cldlayphase(i,ncol,2,5) = 1. 624 else ! low cloud 625 cldlayphase(i,ncol,1,5) = 1. 626 626 endif 627 627 … … 630 630 lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. 631 631 tmpi(i,ncol,nlev)=tmp(i,nlev) 632 633 634 635 636 637 638 632 cldlayphase(i,ncol,4,1) = 1. ! tot cloud 633 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 634 cldlayphase(i,ncol,3,1) = 1. 635 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 636 cldlayphase(i,ncol,2,1) = 1. 637 else ! low cloud 638 cldlayphase(i,ncol,1,1) = 1. 639 639 endif 640 640 … … 651 651 lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. 652 652 tmpl(i,ncol,nlev)=tmp(i,nlev) 653 654 655 cldlayphase(i,ncol,3,2) = 1.656 657 658 659 660 653 cldlayphase(i,ncol,4,2) = 1. ! tot cloud 654 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 655 cldlayphase(i,ncol,3,2) = 1. 656 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 657 cldlayphase(i,ncol,2,2) = 1. 658 else ! low cloud 659 cldlayphase(i,ncol,1,2) = 1. 660 endif 661 661 662 662 else … … 666 666 lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1. ! keep the information "temperature criterium used" 667 667 ! to classify the phase cloud 668 669 670 cldlayphase(i,ncol,3,4) = 1.671 672 673 674 675 676 677 678 cldlayphase(i,ncol,3,1) = 1.679 680 681 682 683 668 cldlayphase(i,ncol,4,4) = 1. ! tot cloud 669 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 670 cldlayphase(i,ncol,3,4) = 1. 671 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 672 cldlayphase(i,ncol,2,4) = 1. 673 else ! low cloud 674 cldlayphase(i,ncol,1,4) = 1. 675 endif 676 cldlayphase(i,ncol,4,1) = 1. ! tot cloud 677 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 678 cldlayphase(i,ncol,3,1) = 1. 679 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 680 cldlayphase(i,ncol,2,1) = 1. 681 else ! low cloud 682 cldlayphase(i,ncol,1,1) = 1. 683 endif 684 684 685 685 endif 686 686 687 687 endif ! end of discrimination condition 688 688 endif ! end of cloud condition 689 689 enddo ! end of altitude loop 690 690 … … 696 696 697 697 toplvlsat=0 698 donlev=17,1,-1 ! from 8.16km until 0km698 DO nlev=17,1,-1 ! from 8.16km until 0km 699 699 p1 = pplay(i,nlev) 700 700 701 701 if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then 702 702 ! Phase discrimination line : ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 703 703 ! + ATB*epsilon50 + zeta50 … … 718 718 lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1. 719 719 720 720 cldlayphase(i,ncol,4,2) = 1. ! tot cloud 721 721 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 722 723 724 725 726 722 cldlayphase(i,ncol,3,2) = 1. 723 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 724 cldlayphase(i,ncol,2,2) = 1. 725 else ! low cloud 726 cldlayphase(i,ncol,1,2) = 1. 727 727 endif 728 728 729 730 731 732 733 734 735 729 cldlayphase(i,ncol,4,5) = 1. ! tot cloud 730 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 731 cldlayphase(i,ncol,3,5) = 1. 732 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 733 cldlayphase(i,ncol,2,5) = 1. 734 else ! low cloud 735 cldlayphase(i,ncol,1,5) = 1. 736 736 endif 737 737 … … 741 741 tmpi(i,ncol,nlev)=tmp(i,nlev) 742 742 743 744 745 746 747 748 749 743 cldlayphase(i,ncol,4,1) = 1. ! tot cloud 744 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 745 cldlayphase(i,ncol,3,1) = 1. 746 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 747 cldlayphase(i,ncol,2,1) = 1. 748 else ! low cloud 749 cldlayphase(i,ncol,1,1) = 1. 750 750 endif 751 751 … … 763 763 tmpl(i,ncol,nlev)=tmp(i,nlev) 764 764 765 766 767 cldlayphase(i,ncol,3,2) = 1.768 769 770 771 772 765 cldlayphase(i,ncol,4,2) = 1. ! tot cloud 766 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 767 cldlayphase(i,ncol,3,2) = 1. 768 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 769 cldlayphase(i,ncol,2,2) = 1. 770 else ! low cloud 771 cldlayphase(i,ncol,1,2) = 1. 772 endif 773 773 774 774 else … … 778 778 lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1. ! false liq ==> ice 779 779 780 781 782 cldlayphase(i,ncol,3,4) = 1.783 784 785 786 787 788 789 790 791 cldlayphase(i,ncol,3,1) = 1.792 793 794 795 796 780 cldlayphase(i,ncol,4,4) = 1. ! tot cloud 781 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 782 cldlayphase(i,ncol,3,4) = 1. 783 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 784 cldlayphase(i,ncol,2,4) = 1. 785 else ! low cloud 786 cldlayphase(i,ncol,1,4) = 1. 787 endif 788 789 cldlayphase(i,ncol,4,1) = 1. ! tot cloud 790 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 791 cldlayphase(i,ncol,3,1) = 1. 792 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 793 cldlayphase(i,ncol,2,1) = 1. 794 else ! low cloud 795 cldlayphase(i,ncol,1,1) = 1. 796 endif 797 797 798 798 endif 799 799 endif ! end of discrimination condition 800 800 801 801 toplvlsat=0 802 802 803 803 ! Find the level of the highest cloud with SR>30 804 if(x(i,ncol,nlev).gt.S_cld_att)then! SR > 30.805 806 goto 99807 808 809 804 if(x(i,ncol,nlev).gt.S_cld_att)then ! SR > 30. 805 toplvlsat=nlev-1 806 goto 99 807 endif 808 809 endif ! end of cloud condition 810 810 enddo ! end of altitude loop 811 811 … … 818 818 !____________________________________________________________________________________________________ 819 819 820 if(toplvlsat.ne.0)then 821 donlev=toplvlsat,1,-1820 if(toplvlsat.ne.0)then 821 DO nlev=toplvlsat,1,-1 822 822 p1 = pplay(i,nlev) 823 823 if(cldy(i,ncol,nlev).eq.1.)then 824 824 lidarcldphase(i,nlev,3)=lidarcldphase(i,nlev,3)+1. 825 825 tmpu(i,ncol,nlev)=tmp(i,nlev) 826 826 827 827 cldlayphase(i,ncol,4,3) = 1. ! tot cloud 828 828 if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud 829 829 cldlayphase(i,ncol,3,3) = 1. 830 830 else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud 831 831 cldlayphase(i,ncol,2,3) = 1. 832 832 else ! low cloud 833 833 cldlayphase(i,ncol,1,3) = 1. 834 835 836 endif 834 endif 835 836 endif 837 837 enddo 838 838 endif … … 876 876 877 877 ! Compute Phase low mid high cloud fractions 878 doiz = 1, Ncat879 doi=1,Nphase-3880 doic = 1, Ncolumns878 DO iz = 1, Ncat 879 DO i=1,Nphase-3 880 DO ic = 1, Ncolumns 881 881 cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i) 882 882 cldlayerphasesum(:,iz)=cldlayerphasesum(:,iz)+cldlayphase(:,ic,iz,i) … … 885 885 enddo 886 886 887 doiz = 1, Ncat888 doi=4,5889 doic = 1, Ncolumns887 DO iz = 1, Ncat 888 DO i=4,5 889 DO ic = 1, Ncolumns 890 890 cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i) 891 891 enddo … … 901 901 ENDWHERE 902 902 903 doi=1,Nphase-1903 DO i=1,Nphase-1 904 904 WHERE ( cldlayerphasesum(:,:).gt.0.0 ) 905 905 cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:) … … 908 908 909 909 910 doi=1,Npoints911 doiz=1,Ncat910 DO i=1,Npoints 911 DO iz=1,Ncat 912 912 checkcldlayerphase=0. 913 913 checkcldlayerphase2=0. 914 914 915 915 if (cldlayerphasesum(i,iz).gt.0.0 )then 916 doic=1,Nphase-3916 DO ic=1,Nphase-3 917 917 checkcldlayerphase=checkcldlayerphase+cldlayerphase(i,iz,ic) 918 918 enddo … … 925 925 enddo 926 926 927 doi=1,Nphase-1927 DO i=1,Nphase-1 928 928 WHERE ( nsublayer(:,:).eq.0.0 ) 929 929 cldlayerphase(:,:,i) = undef … … 934 934 935 935 ! Compute Phase 3D as a function of temperature 936 donlev=1,Nlevels937 do ncol=1,Ncolumns 938 doi=1,Npoints939 doitemp=1,Ntemp936 DO nlev=1,Nlevels 937 DO ncol=1,Ncolumns 938 DO i=1,Npoints 939 DO itemp=1,Ntemp 940 940 if(tmpi(i,ncol,nlev).gt.0.)then 941 941 if( (tmpi(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpi(i,ncol,nlev).lt.tempmod(itemp+1)) )then … … 957 957 958 958 ! Check temperature cloud fraction 959 doi=1,Npoints960 doitemp=1,Ntemp959 DO i=1,Npoints 960 DO itemp=1,Ntemp 961 961 checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4) 962 962 963 964 965 966 963 if(checktemp.NE.lidarcldtemp(i,itemp,1))then 964 print *, i,itemp 965 print *, lidarcldtemp(i,itemp,1:4) 966 endif 967 967 968 968 enddo … … 979 979 ENDWHERE 980 980 981 doi=1,4981 DO i=1,4 982 982 WHERE(lidarcldtempind(:,:).gt.0.) 983 983 lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:) … … 1025 1025 1026 1026 ! #################################################################################### 1027 ! 1) Initialize 1027 ! 1) Initialize 1028 1028 ! #################################################################################### 1029 1029 cldtype = 0.0 … … 1040 1040 ! 2) Cloud detection and Fully attenuated layer detection 1041 1041 ! #################################################################################### 1042 dok=1,Nlevels1042 DO k=1,Nlevels 1043 1043 ! Cloud detection at subgrid-scale: 1044 1044 where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) ) … … 1072 1072 ! #################################################################################### 1073 1073 1074 dok= Nlevels,1,-11075 doic = 1, Ncolumns1076 doip = 1, Npoints1074 DO k= Nlevels,1,-1 1075 DO ic = 1, Ncolumns 1076 DO ip = 1, Npoints 1077 1077 1078 1078 cldlay(ip,ic,1) = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque clouds … … 1090 1090 1091 1091 ! OPAQ variables 1092 doic = 1, Ncolumns1093 doip = 1, Npoints1092 DO ic = 1, Ncolumns 1093 DO ip = 1, Npoints 1094 1094 1095 1095 ! Declaring non-opaque cloudy profiles as thin cloud profiles 1096 1097 1098 1096 if ( (cldlay(ip,ic,4) .eq. 1.0) .and. (cldlay(ip,ic,1) .eq. 0.0) ) then 1097 cldlay(ip,ic,2) = 1.0 1098 endif 1099 1099 1100 1100 ! Filling in 3D and 2D variables 1101 1101 1102 1102 ! Opaque cloud profiles 1103 1104 1105 dok=2,Nlevels1103 if ( cldlay(ip,ic,1) .eq. 1.0 ) then 1104 zopac = 0.0 1105 DO k=2,Nlevels 1106 1106 ! Declaring opaque cloud fraction and z_opaque altitude for 3D and 2D variables 1107 1108 1109 1110 1111 1112 1113 1114 1107 if ( (cldy(ip,ic,k) .eq. 1.0) .and. (zopac .eq. 0.0) ) then 1108 lidarcldtype(ip,k-1,3) = lidarcldtype(ip,k-1,3) + 1.0 1109 cldlay(ip,ic,3) = vgrid_z(k-1) !z_opaque altitude 1110 nsublay(ip,ic,3) = 1.0 1111 zopac = 1.0 1112 endif 1113 if ( cldy(ip,ic,k) .eq. 1.0 ) then 1114 lidarcldtype(ip,k,1) = lidarcldtype(ip,k,1) + 1.0 1115 1115 endif 1116 1117 1116 enddo 1117 endif 1118 1118 1119 1119 ! Thin cloud profiles 1120 1121 dok=1,Nlevels1120 if ( cldlay(ip,ic,2) .eq. 1.0 ) then 1121 DO k=1,Nlevels 1122 1122 ! Declaring thin cloud fraction for 3D variable 1123 1123 if ( cldy(ip,ic,k) .eq. 1.0 ) then 1124 1124 lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1.0 1125 1125 endif 1126 1126 enddo 1127 1127 endif 1128 1128 … … 1146 1146 ! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=Nlevels) to SFC(k=1) 1147 1147 lidarcldtype(:,Nlevels,4) = lidarcldtype(:,Nlevels,3) 1148 doip = 1, Npoints1149 dok = Nlevels-1, 1, -11148 DO ip = 1, Npoints 1149 DO k = Nlevels-1, 1, -1 1150 1150 if ( lidarcldtype(ip,k,3) .ne. undef ) then 1151 1151 lidarcldtype(ip,k,4) = lidarcldtype(ip,k+1,4) + lidarcldtype(ip,k,3) 1152 1152 endif 1153 1153 enddo 1154 1154 enddo 1155 1155 where ( nsubopaq(:,:) .eq. 0.0 ) … … 1159 1159 ! Layered cloud types (opaque, thin and z_opaque 2D variables) 1160 1160 1161 doiz = 1, Ntype1162 doic = 1, Ncolumns1161 DO iz = 1, Ntype 1162 DO ic = 1, Ncolumns 1163 1163 cldtype(:,iz) = cldtype(:,iz) + cldlay(:,ic,iz) 1164 1164 nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_modis_sim.F90
r5099 r5158 200 200 logical, dimension(size(retrievedTau)) :: cloudMask 201 201 real, dimension(size(waterSize, 1), size(waterSize, 2)) :: tauLiquidFraction, tauTotal 202 real:: integratedLiquidFraction202 REAL :: integratedLiquidFraction 203 203 integer :: i, nSubcols, nLevels 204 204 … … 256 256 end if 257 257 258 doi = 1, nSubCols258 DO i = 1, nSubCols 259 259 if(cloudMask(i)) then 260 260 … … 371 371 liquid_opticalThickness, ice_opticalThickness, tauLiquidFraction 372 372 373 real:: seuil373 REAL :: seuil 374 374 ! --------------------------------------------------- 375 375 … … 596 596 reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask) 597 597 reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask) 598 doj=1,nPoints598 DO j=1,nPoints 599 599 600 600 ! Fill clear and optically thin subcolumns with fill … … 651 651 integer :: ij,ik 652 652 653 doij=2,nbin1+1654 doik=2,nbin2+1653 DO ij=2,nbin1+1 654 DO ik=2,nbin2+1 655 655 jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. & 656 656 var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik)) … … 778 778 ! Joint histogram 779 779 780 do i = 1, numTauHistogramBins780 DO i = 1, numTauHistogramBins 781 781 where(cloudMask(:, :)) 782 782 tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. & … … 787 787 end do 788 788 789 do i = 1, numPressureHistogramBins789 DO i = 1, numPressureHistogramBins 790 790 where(cloudMask(:, :)) 791 791 pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. & … … 796 796 end do 797 797 798 doi = 1, numPressureHistogramBins799 doj = 1, numTauHistogramBins798 DO i = 1, numPressureHistogramBins 799 DO j = 1, numTauHistogramBins 800 800 Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = & 801 801 real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols) … … 808 808 real, dimension(:), intent(in) :: tauIncrement, pressure 809 809 real, intent(in) :: tauLimit 810 real:: cloud_top_pressure810 REAL :: cloud_top_pressure 811 811 812 812 ! Find the extinction-weighted pressure. Assume that pressure varies linearly between 813 813 ! layers and use the trapezoidal rule. 814 814 815 real:: deltaX, totalTau, totalProduct815 REAL :: deltaX, totalTau, totalProduct 816 816 integer :: i 817 817 818 818 totalTau = 0.; totalProduct = 0. 819 doi = 2, size(tauIncrement)819 DO i = 2, size(tauIncrement) 820 820 if(totalTau + tauIncrement(i) > tauLimit) then 821 821 deltaX = tauLimit - totalTau … … 840 840 real, dimension(:), intent(in) :: tauIncrement, f 841 841 real, intent(in) :: tauLimit 842 real:: weight_by_extinction842 REAL :: weight_by_extinction 843 843 844 844 ! Find the extinction-weighted value of f(tau), assuming constant f within each layer 845 845 846 real:: deltaX, totalTau, totalProduct846 REAL :: deltaX, totalTau, totalProduct 847 847 integer :: i 848 848 849 849 totalTau = 0.; totalProduct = 0. 850 doi = 1, size(tauIncrement)850 DO i = 1, size(tauIncrement) 851 851 if(totalTau + tauIncrement(i) > tauLimit) then 852 852 deltaX = tauLimit - totalTau … … 864 864 pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size) 865 865 real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size 866 real:: compute_nir_reflectance866 REAL :: compute_nir_reflectance 867 867 868 868 real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, & … … 893 893 integer, intent(in) :: phase 894 894 real, intent(in) :: tau, obs_Refl_nir 895 real:: retrieve_re895 REAL :: retrieve_re 896 896 897 897 ! Finds the re that produces the minimum mis-match between predicted and observed reflectance in … … 907 907 908 908 real, parameter :: min_distance_to_boundary = 0.01 909 real:: re_min, re_max, delta_re909 REAL :: re_min, re_max, delta_re 910 910 integer :: i 911 911 … … 957 957 real, dimension(:), intent(in) :: x, y 958 958 real, intent(in) :: yobs 959 real:: interpolate_to_min959 REAL :: interpolate_to_min 960 960 961 961 ! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs) … … 1008 1008 integer, intent(in) :: phase 1009 1009 real, intent(in) :: re 1010 real :: get_g_nir1010 REAL :: get_g_nir 1011 1011 1012 1012 real, dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), & … … 1035 1035 integer, intent(in) :: phase 1036 1036 real, intent(in) :: re 1037 real:: get_ssa_nir1037 REAL :: get_ssa_nir 1038 1038 1039 1039 ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function … … 1060 1060 real, intent(in) :: x 1061 1061 real, dimension(:), intent(in) :: coefficients 1062 real:: fit_to_cubic1062 REAL :: fit_to_cubic 1063 1063 1064 1064 … … 1069 1069 real, intent(in) :: x 1070 1070 real, dimension(:), intent(in) :: coefficients 1071 real:: fit_to_quadratic1071 REAL :: fit_to_quadratic 1072 1072 1073 1073 … … 1079 1079 pure function compute_toa_reflectace(tau, g, w0) 1080 1080 real, dimension(:), intent(in) :: tau, g, w0 1081 real:: compute_toa_reflectace1081 REAL :: compute_toa_reflectace 1082 1082 1083 1083 logical, dimension(size(tau)) :: cloudMask 1084 1084 integer, dimension(count(tau(:) > 0)) :: cloudIndicies 1085 1085 real, dimension(count(tau(:) > 0)) :: Refl, Trans 1086 real:: Refl_tot, Trans_tot1086 REAL :: Refl_tot, Trans_tot 1087 1087 integer :: i 1088 1088 ! --------------------------------------- … … 1092 1092 cloudMask = tau(:) > 0. 1093 1093 cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask) 1094 doi = 1, size(cloudIndicies)1094 DO i = 1, size(cloudIndicies) 1095 1095 call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i)) 1096 1096 end do … … 1115 1115 integer, parameter :: beam = 2 1116 1116 real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999 1117 real:: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &1117 REAL :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & 1118 1118 rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th 1119 1119 … … 1183 1183 elemental function two_stream_reflectance(tauint, gint, w0int) 1184 1184 real, intent(in) :: tauint, gint, w0int 1185 real:: two_stream_reflectance1185 REAL :: two_stream_reflectance 1186 1186 1187 1187 ! Compute reflectance in a single layer using the two stream approximation … … 1194 1194 integer, parameter :: beam = 2 1195 1195 real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999 1196 real:: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &1196 REAL :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & 1197 1197 rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den 1198 1198 ! ------------------------ … … 1267 1267 Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1) 1268 1268 1269 doi=2, size(Refl)1269 DO i=2, size(Refl) 1270 1270 ! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface): 1271 1271 Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1 - Refl_cumulative(i-1) * Refl(i)) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/optics_lib.F90
r5099 r5158 519 519 520 520 ! // region from 0.045 microns to 167.0 microns - no temperature depend 521 doi=2,nwl521 DO i=2,nwl 522 522 if(alam < wl(i)) continue 523 523 enddo … … 539 539 if(tk > temref(1)) tk=temref(1) 540 540 if(tk < temref(4)) tk=temref(4) 541 do11 i=2,4541 DO 11 i=2,4 542 542 if(tk.ge.temref(i)) go to 12 543 543 11 continue 544 544 12 lt1=i 545 545 lt2=i-1 546 do13 i=2,nwlt546 DO 13 i=2,nwlt 547 547 if(alam.le.wlt(i)) go to 14 548 548 13 continue -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/phys_cosp.F90
r5133 r5158 44 44 ! mr_ozone, !Concentration ozone (Kg/Kg) 45 45 ! dem_s !Cloud optical emissivity 46 ! dtau_s 47 ! emsfc_lw = 1. 46 ! dtau_s !Cloud optical thickness 47 ! emsfc_lw = 1. !Surface emissivity dans radlwsw.F90 48 48 49 49 !!! Outputs : … … 132 132 ! Declaration necessaires pour les sorties IOIPSL 133 133 integer :: ii 134 real:: ecrit_day,ecrit_hf,ecrit_mth, missing_val134 REAL :: ecrit_day,ecrit_hf,ecrit_mth, missing_val 135 135 logical :: ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, ok_all_xml 136 136 … … 150 150 real,dimension(Nlevlmdz) :: presnivs 151 151 integer :: itap,k,ip 152 real:: dtime,freq_cosp152 REAL :: dtime,freq_cosp 153 153 real,dimension(2) :: time_bnds 154 154 … … 249 249 250 250 zlev_half(:,1) = phis(:)/9.81 251 dok = 2, Nlevels252 doip = 1, Npoints251 DO k = 2, Nlevels 252 DO ip = 1, Npoints 253 253 zlev_half(ip,k) = phi(ip,k)/9.81 + & 254 254 (phi(ip,k)-phi(ip,k-1))/9.81 * (ph(ip,k)-p(ip,k)) / (p(ip,k)-p(ip,k-1)) … … 267 267 gbx%skt = skt !Skin temperature (K) 268 268 269 doip = 1, Npoints269 DO ip = 1, Npoints 270 270 if (fracTerLic(ip).ge.0.5) then 271 271 gbx%land(ip) = 1. -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator.F90
r5099 r5158 41 41 ! [hm_matrix] table of hydrometeor mixing rations (g/kg) 42 42 ! [re_matrix] table of hydrometeor effective radii. 0 ==> use defaults. (units=microns) 43 ! [Np_matrix] 43 ! [Np_matrix] table of hydrometeor number concentration. 0 ==> use defaults. (units = 1/kg) 44 44 45 45 ! Outputs: … … 117 117 ns ! number of discrete drop sizes 118 118 119 logical :: hydro 119 logical :: hydro ! true=hydrometeor in vol, false=none 120 120 real*8 :: & 121 121 rho_a, & ! air density (kg m^-3) … … 180 180 d_gate=-1 181 181 endif 182 dok=start_gate,end_gate,d_gate182 DO k=start_gate,end_gate,d_gate 183 183 ! // loop over each profile (nprof) 184 dopr=1,nprof184 DO pr=1,nprof 185 185 t_kelvin = t_matrix(pr,k) 186 186 ! :: determine if hydrometeor(s) present in volume 187 187 hydro = .false. 188 doj=1,hp%nhclass188 DO j=1,hp%nhclass 189 189 if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then 190 190 hydro = .true. … … 198 198 rho_a = (p_matrix(pr,k)*100.)/(287.0*(t_kelvin)) 199 199 ! :: loop over hydrometeor type 200 dotp=1,hp%nhclass200 DO tp=1,hp%nhclass 201 201 if (hm_matrix(tp,pr,k) <= 1E-12) cycle 202 202 phase = hp%phase(tp) … … 233 233 if (iRe_type.lt.1) iRe_type=1 234 234 235 Re=step*(iRe_type+0.5) 235 Re=step*(iRe_type+0.5) ! set value of Re to closest value allowed in LUT. 236 236 iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step) 237 237 … … 273 273 if (hp%rho(tp) < 0) then 274 274 ! Use equivalent volume spheres. 275 hp%rho_eff(tp,1:ns,iRe_type) = 917 275 hp%rho_eff(tp,1:ns,iRe_type) = 917 ! solid ice == equivalent volume approach 276 276 Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * ( (Di*1E-6) ** (hp%bpm(tp)/3.0) ) * 1E6 277 277 ! alternative is to comment out above two lines and use the following block … … 362 362 endif 363 363 364 enddo 364 enddo ! end loop of tp (hydrometeor type) 365 365 366 366 else 367 ! :: volume is hydrometeor-free 367 ! :: volume is hydrometeor-free 368 368 kr_vol(pr,k) = 0 369 369 z_vol(pr,k) = undef -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator_init.F90
r5099 r5158 60 60 integer, intent(in) :: nhclass ! number of hydrometeor classes in use 61 61 integer, intent(in) :: use_gas_abs,do_ray 62 real:: undef62 REAL :: undef 63 63 real,dimension(nhclass),intent(in) :: hclass_dmin,hclass_dmax, & 64 64 hclass_apm,hclass_bpm,hclass_rho, & … … 90 90 ! Store settings for hydrometeor types in hp (class_parameter) structure. 91 91 92 doi = 1,nhclass92 DO i = 1,nhclass 93 93 hp%dtype(i) = hclass_type(i) 94 94 hp%phase(i) = hclass_phase(i) … … 117 117 118 118 hp%base_list(1)=0; 119 doj=1,Re_MAX_BIN119 DO j=1,Re_MAX_BIN 120 120 hp%step_list(j)=0.1+0.1*((j-1)**1.5); 121 121 if(hp%step_list(j)>Re_BIN_LENGTH) then … … 129 129 ! set up Temperature bin structure used for z scaling 130 130 131 doi=1,cnt_ice131 DO i=1,cnt_ice 132 132 hp%mt_tti(i)=(i-1)*5-90 + 273.15 133 133 enddo 134 134 135 doi=1,cnt_liq135 DO i=1,cnt_liq 136 136 hp%mt_ttl(i)=(i-1)*5-60 + 273.15 137 137 enddo … … 143 143 144 144 hp%D(1) = dmin 145 doi=2,nd145 DO i=2,nd 146 146 hp%D(i) = hp%D(i-1)*deltp 147 147 enddo -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/radar_simulator_types.F90
r5095 r5158 45 45 46 46 ! defines location of radar relative to hgt_matrix. 47 logical :: radar_at_layer_one ! if true radar is assume to be at the edge48 49 50 51 ! a space borne radar.47 logical :: radar_at_layer_one ! if true radar is assume to be at the edge 48 ! of the first layer, if the first layer is the 49 ! surface than a ground-based radar. If the 50 ! first layer is the top-of-atmosphere, then 51 ! a space borne radar. 52 52 53 53 ! variables used to store Z scale factors -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/scale_luts_io.F90
r5099 r5158 42 42 trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' 43 43 44 doi=1,maxhclass45 doj=1,mt_ntt46 dok=1,nRe_types44 DO i=1,maxhclass 45 DO j=1,mt_ntt 46 DO k=1,nRe_types 47 47 48 48 ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt) … … 95 95 trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' 96 96 97 doi=1,maxhclass98 doj=1,mt_ntt99 dok=1,nRe_types97 DO i=1,maxhclass 98 DO j=1,mt_ntt 99 DO k=1,nRe_types 100 100 101 101 ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/zeff.F90
r5099 r5158 110 110 sizep = (pi*D0)/wl 111 111 dqv(1) = 0. 112 doi=1,nsizes112 DO i=1,nsizes 113 113 call mieint(sizep(i), m0(i), one, dqv, qext(i), dqsc, qbsca(i), & 114 114 dg, xs1, xs2, dph, err)
Note: See TracChangeset
for help on using the changeset viewer.