Changeset 5081 for LMDZ6/branches/Amaury_dev/libf/phylmd
- Timestamp:
- Jul 19, 2024, 4:15:44 PM (6 months ago)
- Location:
- LMDZ6/branches/Amaury_dev/libf/phylmd
- Files:
-
- 30 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/cond_evap_tstep_mod.F90
r4950 r5081 227 227 ! H2SO4 mass fraction in aerosol 228 228 WH2=R2SO4*1.0E-2 229 IF(WH2 .EQ.0.0) RETURN229 IF(WH2==0.0) RETURN 230 230 ! ACTIVITY COEFFICIENT(SEE GIAUQUE,1951) 231 231 ! AYERS ET AL (1980) … … 324 324 ff(IK,:)=0.0 325 325 DO k=1, nbtr_bin 326 IF (k .LE.(nbtr_bin-1)) THEN327 IF (Vbin_wet(k) .LE.Vnew.AND.Vnew.LT.Vbin_wet(k+1)) THEN326 IF (k<=(nbtr_bin-1)) THEN 327 IF (Vbin_wet(k)<=Vnew.AND.Vnew<Vbin_wet(k+1)) THEN 328 328 ff(IK,k)= Vbin_wet(k)/Vnew*(Vbin_wet(k+1)-Vnew)/(Vbin_wet(k+1)-Vbin_wet(k)) 329 329 ENDIF 330 330 ENDIF 331 IF (k .EQ.1.AND.Vnew.LE.Vbin_wet(k)) THEN331 IF (k==1.AND.Vnew<=Vbin_wet(k)) THEN 332 332 ff(IK,k)= 1. 333 333 ENDIF 334 IF (k .GT.1) THEN335 IF (Vbin_wet(k-1) .LT.Vnew.AND.Vnew.LT.Vbin_wet(k)) THEN334 IF (k>1) THEN 335 IF (Vbin_wet(k-1)<Vnew.AND.Vnew<Vbin_wet(k)) THEN 336 336 ff(IK,k)= 1.-ff(IK,k-1) 337 337 ENDIF 338 338 ENDIF 339 IF (k .EQ.nbtr_bin.AND.Vnew.GE.Vbin_wet(k)) THEN339 IF (k==nbtr_bin.AND.Vnew>=Vbin_wet(k)) THEN 340 340 ff(IK,k)= 1. 341 341 ENDIF -
LMDZ6/branches/Amaury_dev/libf/phylmd/YOECUMF.h
r1907 r5081 12 12 ! 13 13 COMMON /YOECUMF/ & 14 & LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV, &15 14 & ENTRPEN,ENTRSCV,ENTRMID,ENTRDD,CMFCTOP, & 16 & CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON 15 & CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON, & 16 & LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV 17 17 18 18 19 LOGICAL LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV -
LMDZ6/branches/Amaury_dev/libf/phylmd/aaam_bud.F90
r2350 r5081 127 127 hadley = 1.E18 128 128 hadday = 1.E18*24.*3600. 129 IF(klon_glo .EQ.1) THEN129 IF(klon_glo==1) THEN 130 130 dlat = xpi 131 131 ELSE -
LMDZ6/branches/Amaury_dev/libf/phylmd/acama_gwd_rando_m.F90
r3977 r5081 358 358 CMAX*2.*(MOD(TT(II, JW+4*(JJ-1)+JJ)**2, 1.)-0.5)*SQRT(3.)/SQRT(NA*1.) 359 359 END DO 360 IF (CPHA .LT.0.) THEN360 IF (CPHA<0.) THEN 361 361 CPHA = -1.*CPHA 362 362 ZP(JW,II) = ZP(JW,II) + RPI -
LMDZ6/branches/Amaury_dev/libf/phylmd/add_phys_tend_mod.F90
r5051 r5081 226 226 ENDDO 227 227 228 IF (fl_ebil .GT.0) THEN228 IF (fl_ebil > 0) THEN 229 229 ! ------------------------------------------------ 230 230 ! Compute vertical sum for each atmospheric column … … 283 283 !===================================================================================== 284 284 285 IF (jbad .GT.0) THEN285 IF (jbad > 0) THEN 286 286 DO j = 1, jbad 287 287 i=jadrs(j) 288 IF (prt_level .ge.debug_level) THEN288 IF (prt_level>=debug_level) THEN 289 289 print*,'PLANTAGE POUR LE POINT i lon lat =',& 290 290 i,longitude_deg(i),latitude_deg(i),text … … 301 301 ! Impression, warning et correction en cas de probleme moins important 302 302 !===================================================================================== 303 IF (jqbad .GT.0) THEN303 IF (jqbad > 0) THEN 304 304 done(:) = .false. !jyg 305 305 DO j = 1, jqbad 306 306 i=jqadrs(j) 307 if(prt_level .ge.debug_level) THEN307 if(prt_level>=debug_level) THEN 308 308 print*,'WARNING : EAU POUR LE POINT i lon lat =',& 309 309 i,longitude_deg(i),latitude_deg(i),text … … 325 325 zqp_int = zqp_int + zqp(k) *(paprs(i,k)-paprs(i,k+1))/Rg 326 326 ENDDO 327 IF (prt_level .ge.debug_level) THEN327 IF (prt_level>=debug_level) THEN 328 328 print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', & 329 329 i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int … … 340 340 DO k = 1, klev 341 341 zq=q_seri(i,k)+zdq(i,k) 342 IF (zq .lt.1.e-15) THEN343 IF (q_seri(i,k) .lt.1.e-15) THEN344 IF (prt_level .ge.debug_level) THEN342 IF (zq<1.e-15) THEN 343 IF (q_seri(i,k)<1.e-15) THEN 344 IF (prt_level>=debug_level) THEN 345 345 print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k) 346 346 ENDIF … … 383 383 ENDDO 384 384 ENDDO 385 IF (jbad .GT.0) THEN385 IF (jbad > 0) THEN 386 386 DO j = 1, jbad 387 387 i=jadrs(j) 388 388 k=kadrs(j) 389 if(prt_level .ge.debug_level) THEN389 if(prt_level>=debug_level) THEN 390 390 print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',& 391 391 i,itap,longitude_deg(i),latitude_deg(i),text,jbad, & … … 401 401 ENDIF 402 402 ! 403 IF (jqbad .GT.0) THEN403 IF (jqbad > 0) THEN 404 404 DO j = 1, jqbad 405 405 i=jqadrs(j) 406 406 k=kqadrs(j) 407 IF (prt_level .ge.debug_level) THEN407 IF (prt_level>=debug_level) THEN 408 408 print*,'WARNING : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',& 409 409 i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,& … … 441 441 !====================================================================== 442 442 443 IF (fl_ebil .GT.0) THEn443 IF (fl_ebil > 0) THEn 444 444 445 445 ! ------------------------------------------------ … … 584 584 ENDDO 585 585 586 IF (fl_ebil .GT.0) THEN586 IF (fl_ebil > 0) THEN 587 587 ! ------------------------------------------------ 588 588 ! Compute vertical sum for each atmospheric column … … 613 613 !====================================================================== 614 614 615 IF (fl_ebil .GT.0) THEN615 IF (fl_ebil > 0) THEN 616 616 617 617 ! ------------------------------------------------ … … 761 761 762 762 !!print *,'prt_level:',prt_level,' fl_ebil:',fl_ebil,' fl_cor_ebil:',fl_cor_ebil 763 IF ((fl_ebil .GT. 0) .AND. (klon .EQ.1)) THEN763 IF ((fl_ebil > 0) .AND. (klon == 1)) THEN 764 764 765 765 bilq_bnd = 0. … … 793 793 bilh_error = d_h_col(1) - bilh_bnd 794 794 ! are the errors too large? 795 IF (abs(bilq_error) .GT.bilq_seuil) bilq_ok=1796 IF (abs(bilh_error) .GT.bilh_seuil) bilh_ok=1795 IF (abs(bilq_error) > bilq_seuil) bilq_ok=1 796 IF (abs(bilh_error) > bilh_seuil) bilh_ok=1 797 797 ! 798 798 ! Print diagnostics 799 799 ! ================= 800 IF ( (bilq_ok .eq. 0).AND.(bilh_ok .eq.0) ) THEN800 IF ( (bilq_ok == 0).AND.(bilh_ok == 0) ) THEN 801 801 status="enerbil-OK" 802 802 ELSE … … 804 804 ENDIF 805 805 806 IF (prt_level .GE.3) THEN806 IF (prt_level >= 3) THEN 807 807 write(*,9010) text,status," itap:",itap,"enerbilERROR: Q", bilq_error," H", bilh_error 808 808 9010 format (1x,A8,2x,A12,A6,I4,A18,E15.6,A5,E15.6) 809 809 ENDIF 810 IF (prt_level .GE.3) THEN810 IF (prt_level >= 3) THEN 811 811 write(*,9000) text,"enerbil: Q,H,KE budget", d_qt_col(1), d_h_col(1),d_ek_col(1) 812 812 ENDIF 813 IF (prt_level .GE.5) THEN813 IF (prt_level >= 5) THEN 814 814 write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd 815 815 write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1), d_qbs_col(1) … … 819 819 specific_diag: SELECT CASE (text) 820 820 CASE("vdf") specific_diag 821 IF (prt_level .GE.5) THEN821 IF (prt_level >= 5) THEN 822 822 write(*,9000) text,"enerbil: d_h, bilh, sens,t_seri", d_h_col(1), bilh_bnd, sens(1), t_seri(1,1) 823 823 write(*,9000) text,"enerbil: d_h_col_vdf, f_h, diff",d_h_col_vdf, f_h_bnd, bilh_bnd-sens(1) 824 824 ENDIF 825 825 CASE("lsc") specific_diag 826 IF (prt_level .GE.5) THEN826 IF (prt_level >= 5) THEN 827 827 write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_lsc(1), rlvtt * rain_lsc(1), -(rcw-rcpd)*t_seri(1,1) * rain_lsc(1) 828 828 write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_lsc(1), rlstt * snow_lsc(1), -(rcs-rcpd)*t_seri(1,1) * snow_lsc(1) 829 829 ENDIF 830 830 CASE("convection") specific_diag 831 IF (prt_level .GE.5) THEN831 IF (prt_level >= 5) THEN 832 832 write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_con(1), rlvtt * rain_con(1), -(rcw-rcpd)*t_seri(1,1) * rain_con(1) 833 833 write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_con(1), rlstt * snow_con(1), -(rcs-rcpd)*t_seri(1,1) * snow_con(1) -
LMDZ6/branches/Amaury_dev/libf/phylmd/add_wake_tend.F90
r4744 r5081 1 1 SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zdas, zddensw, zddensaw, zoccur, text, abortphy) 2 2 !=================================================================== 3 ! Ajoute les tendances li ées aux diverses parametrisations physiques aux3 ! Ajoute les tendances li�es aux diverses parametrisations physiques aux 4 4 ! variables d'etat des poches froides. 5 5 !=================================================================== … … 43 43 DO l = 1, klev 44 44 DO i = 1, klon 45 IF (zoccur(i) .GE.1) THEN45 IF (zoccur(i) >= 1) THEN 46 46 wake_deltat(i, l) = wake_deltat(i, l) + zddeltat(i,l) 47 47 wake_deltaq(i, l) = wake_deltaq(i, l) + zddeltaq(i,l) … … 53 53 END DO 54 54 DO i = 1, klon 55 IF (zoccur(i) .GE.1) THEN55 IF (zoccur(i) >= 1) THEN 56 56 wake_s(i) = wake_s(i) + zds(i) 57 57 awake_s(i) = awake_s(i) + zdas(i) -
LMDZ6/branches/Amaury_dev/libf/phylmd/alpale.h
r4843 r5081 14 14 15 15 common/calpale1/iflag_trig_bl,iflag_clos_bl,tau_trig_shallow,tau_trig_deep,iflag_strig 16 common/calpale2/ s_trig,iflag_coupl,iflag_clos,iflag_wake,alp_bl_k,h_trig16 common/calpale2/alp_bl_k,s_trig,h_trig,iflag_coupl,iflag_clos,iflag_wake 17 17 18 18 !$OMP THREADPRIVATE(/calpale1/,/calpale2/) -
LMDZ6/branches/Amaury_dev/libf/phylmd/alpale_wk.F90
r3068 r5081 58 58 !! print *,'alpale_wk: sigmaw(1), wdens(1) ', sigmaw(1), wdens(1) 59 59 DO i = 1,klon 60 IF (zoccur(i) .GE.1) THEN60 IF (zoccur(i) >= 1) THEN 61 61 wkrad(i) = sqrt(sigmaw(i)/(rpi*wdens(i))) 62 62 ELSE -
LMDZ6/branches/Amaury_dev/libf/phylmd/clesphys.h
r5017 r5081 133 133 & , qsol0,albsno0,evap0 & 134 134 & , co2_ppm0 & 135 & , tau_thermals & 135 136 !FC 136 137 & , Cd_frein,zrel_oro_t,zpmm_orodr_t,zpmm_orolf_t,zstd_orodr_t & … … 163 164 & , ok_daily_climoz, ok_all_xml, ok_lwoff & 164 165 & , iflag_phytrac, ok_new_lscp, ok_bs, ok_rad_bs & 165 & , iflag_thermals,nsplit_thermals , tau_thermals&166 & , iflag_thermals,nsplit_thermals & 166 167 & , iflag_physiq, ok_3Deffect, ok_water_mass_fixer 167 168 save /clesphys/ -
LMDZ6/branches/Amaury_dev/libf/phylmd/climb_hq_mod.F90
r3102 r5081 7 7 8 8 IMPLICIT NONE 9 SAVE10 9 PRIVATE 11 10 PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd … … 222 221 !!! jyg le 07/02/2012 223 222 !!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN 224 IF (mod(iflag_pbl_split,10) .ge.1) THEN223 IF (mod(iflag_pbl_split,10) >=1) THEN 225 224 !!! nrlmd le 02/05/2011 226 225 DO k= 1, klev … … 231 230 Dcoef_Q_out(i,k) = Dcoef_Q(i,k) 232 231 Kcoef_hq_out(i,k) = Kcoefhq(i,k) 233 IF (k .eq.1) THEN232 IF (k==1) THEN 234 233 gama_h_out(i,k) = 0. 235 234 gama_q_out(i,k) = 0. … … 379 378 !!! jyg le 07/02/2012 380 379 !!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN 381 IF (mod(iflag_pbl_split,10) .ge.1) THEN380 IF (mod(iflag_pbl_split,10) >=1) THEN 382 381 !!! nrlmd le 02/05/2011 383 382 DO i = 1, knon … … 394 393 Dcoef_Q(i,k)=Dcoef_Q_in(i,k) 395 394 Kcoefhq(i,k)=Kcoef_hq_in(i,k) 396 IF (k .gt.1) THEN395 IF (k>1) THEN 397 396 gamah(i,k)=gama_h_in(i,k) 398 397 gamaq(i,k)=gama_q_in(i,k) -
LMDZ6/branches/Amaury_dev/libf/phylmd/conf_phys_m.F90
r4843 r5081 1352 1352 !Config Help = 1353 1353 ! 1354 ok_ice_sursat_omp = 01354 ok_ice_sursat_omp = .FALSE. 1355 1355 CALL getin('ok_ice_sursat',ok_ice_sursat_omp) 1356 1356 … … 2574 2574 2575 2575 !--test on radiative scheme 2576 IF (iflag_rrtm .EQ. 0) THEN2577 IF (NSW .NE.2) THEN2576 IF (iflag_rrtm == 0) THEN 2577 IF (NSW/=2) THEN 2578 2578 WRITE(lunout,*) ' ERROR iflag_rrtm=0 and NSW<>2 not possible' 2579 2579 CALL abort_physic('conf_phys','choice NSW not valid',1) 2580 2580 ENDIF 2581 ELSE IF (iflag_rrtm .EQ.1) THEN2582 IF (NSW .NE.2.AND.NSW.NE.4.AND.NSW.NE.6) THEN2581 ELSE IF (iflag_rrtm == 1) THEN 2582 IF (NSW/=2.AND.NSW/=4.AND.NSW/=6) THEN 2583 2583 WRITE(lunout,*) ' ERROR iflag_rrtm=1 and NSW<>2,4,6 not possible' 2584 2584 CALL abort_physic('conf_phys','choice NSW not valid',1) 2585 2585 ENDIF 2586 ELSE IF (iflag_rrtm .EQ.2) THEN2587 IF (NSW .NE.2.AND.NSW.NE.4.AND.NSW.NE.6) THEN2586 ELSE IF (iflag_rrtm == 2) THEN 2587 IF (NSW/=2.AND.NSW/=4.AND.NSW/=6) THEN 2588 2588 WRITE(lunout,*) ' ERROR iflag_rrtm=1 and NSW<>2,4,6 not possible' 2589 2589 CALL abort_physic('conf_phys','choice NSW not valid',1) … … 2610 2610 2611 2611 !--test on ocean surface albedo 2612 IF (iflag_albedo .LT.0.OR.iflag_albedo.GT.2) THEN2612 IF (iflag_albedo<0.OR.iflag_albedo>2) THEN 2613 2613 WRITE(lunout,*) ' ERROR iflag_albedo<>0,1' 2614 2614 CALL abort_physic('conf_phys','choice iflag_albedo not valid',1) … … 2617 2617 ! Flag_aerosol cannot be set to zero if aerosol direct effect (ade) or aerosol indirect effect (aie) are activated 2618 2618 IF (ok_ade .OR. ok_aie) THEN 2619 IF ( flag_aerosol .EQ.0 ) THEN2619 IF ( flag_aerosol == 0 ) THEN 2620 2620 CALL abort_physic('conf_phys','flag_aerosol=0 not compatible avec ok_ade ou ok_aie=.TRUE.',1) 2621 2621 ENDIF … … 2623 2623 2624 2624 ! Flag_aerosol cannot be set to zero if we are in coupled mode for aerosol 2625 IF (aerosol_couple .AND. flag_aerosol .EQ.0 ) THEN2625 IF (aerosol_couple .AND. flag_aerosol == 0 ) THEN 2626 2626 CALL abort_physic('conf_phys', 'flag_aerosol cannot be to zero if aerosol_couple=y ', 1) 2627 2627 ENDIF 2628 2628 2629 2629 ! Read_climoz needs to be set zero if we are in couple mode for chemistry 2630 IF (chemistry_couple .AND. read_climoz .ne. 0) THEN2630 IF (chemistry_couple .AND. read_climoz /= 0) THEN 2631 2631 CALL abort_physic('conf_phys', 'read_climoz need to be to zero if chemistry_couple=y ', 1) 2632 2632 ENDIF 2633 2633 2634 2634 ! flag_aerosol need to be different to zero if ok_cdnc is activated 2635 IF (ok_cdnc .AND. flag_aerosol .EQ.0) THEN2635 IF (ok_cdnc .AND. flag_aerosol == 0) THEN 2636 2636 CALL abort_physic('conf_phys', 'flag_aerosol cannot be to zero if ok_cdnc is activated ', 1) 2637 2637 ENDIF … … 2643 2643 2644 2644 ! flag_aerosol=7 => MACv2SP climatology 2645 IF (flag_aerosol .EQ.7.AND. iflag_rrtm.NE.1) THEN2645 IF (flag_aerosol==7.AND. iflag_rrtm/=1) THEN 2646 2646 CALL abort_physic('conf_phys', 'flag_aerosol=7 (MACv2SP) can only be activated with RRTM',1) 2647 2647 ENDIF 2648 IF (flag_aerosol .EQ.7.AND. NSW.NE.6) THEN2648 IF (flag_aerosol==7.AND. NSW/=6) THEN 2649 2649 CALL abort_physic('conf_phys', 'flag_aerosol=7 (MACv2SP) can only be activated with NSW=6',1) 2650 2650 ENDIF 2651 2651 2652 2652 ! BC internal mixture is only possible with RRTM & NSW=6 & flag_aerosol=6 or aerosol_couple 2653 IF (flag_bc_internal_mixture .AND. NSW .NE.6) THEN2653 IF (flag_bc_internal_mixture .AND. NSW/=6) THEN 2654 2654 CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with NSW=6',1) 2655 2655 ENDIF 2656 IF (flag_bc_internal_mixture .AND. iflag_rrtm .NE.1) THEN2656 IF (flag_bc_internal_mixture .AND. iflag_rrtm/=1) THEN 2657 2657 CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with RRTM',1) 2658 2658 ENDIF 2659 IF (flag_bc_internal_mixture .AND. flag_aerosol .NE.6) THEN2659 IF (flag_bc_internal_mixture .AND. flag_aerosol/=6) THEN 2660 2660 CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with flag_aerosol=6',1) 2661 2661 ENDIF 2662 2662 2663 2663 ! test sur flag_volc_surfstrat 2664 IF (flag_volc_surfstrat .LT.0.OR.flag_volc_surfstrat.GT.2) THEN2664 IF (flag_volc_surfstrat<0.OR.flag_volc_surfstrat>2) THEN 2665 2665 CALL abort_physic('conf_phys', 'flag_volc_surfstrat can only be 0 1 or 2',1) 2666 2666 ENDIF 2667 IF ((.NOT.ok_volcan.OR..NOT.ok_ade.OR..NOT.ok_aie).AND.flag_volc_surfstrat .GT.0) THEN2667 IF ((.NOT.ok_volcan.OR..NOT.ok_ade.OR..NOT.ok_aie).AND.flag_volc_surfstrat>0) THEN 2668 2668 CALL abort_physic('conf_phys', 'ok_ade, ok_aie, ok_volcan need to be activated if flag_volc_surfstrat is 1 or 2',1) 2669 2669 ENDIF … … 2679 2679 2680 2680 ! Test on chemistry cycle 2681 IF ((type_trac .ne. "inca" .AND. type_trac .ne. "inco") .AND. ( dms_cycle_cpl .OR. n2o_cycle_cpl) ) THEN2681 IF ((type_trac /= "inca" .AND. type_trac /= "inco") .AND. ( dms_cycle_cpl .OR. n2o_cycle_cpl) ) THEN 2682 2682 CALL abort_physic('conf_phys', 'dms_cycle_cpl or n2o_cycle_cpl has to be TRUE only with INCA coupling model',1) 2683 2683 ENDIF -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/m_mrgrnk.F90
r3241 r5081 33 33 IRNGT (1) = 1 34 34 Return 35 Case Default36 Continue37 35 End Select 38 36 ! … … 233 231 IRNGT (1) = 1 234 232 Return 235 Case Default236 Continue237 233 End Select 238 234 ! -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/optics_lib.F90
r2428 r5081 38 38 39 39 ! ----- INPUTS ----- 40 real *8, intent(in) :: freq,tk40 real(kind=8), intent(in) :: freq,tk 41 41 42 42 ! ----- OUTPUTS ----- 43 real *8, intent(out) :: n_r, n_i43 real(kind=8), intent(out) :: n_r, n_i 44 44 45 45 ! ----- INTERNAL ----- 46 real *8ld,es,ei,a,ls,sg,tm1,cos1,sin147 real *8e_r,e_i48 real *8pi49 real *8tc50 complex *16e_comp, sq46 real(kind=8) ld,es,ei,a,ls,sg,tm1,cos1,sin1 47 real(kind=8) e_r,e_i 48 real(kind=8) pi 49 real(kind=8) tc 50 complex(kind=8) e_comp, sq 51 51 52 52 tc = tk - 273.15 … … 102 102 103 103 ! ----- INPUTS ----- 104 real *8, intent(in) :: freq, t104 real(kind=8), intent(in) :: freq, t 105 105 106 106 ! ----- OUTPUTS ----- 107 real *8, intent(out) :: n_r,n_i107 real(kind=8), intent(out) :: n_r,n_i 108 108 109 109 ! Parameters: 110 integer *2:: i,lt1,lt2,nwl,nwlt110 integer(kind=2) :: i,lt1,lt2,nwl,nwlt 111 111 parameter(nwl=468,nwlt=62) 112 112 113 real *8:: alam,cutice,pi,t1,t2,wlmax,wlmin, &113 real(kind=8) :: alam,cutice,pi,t1,t2,wlmax,wlmin, & 114 114 x,x1,x2,y,y1,y2,ylo,yhi,tk 115 115 116 real *8:: &116 real(kind=8) :: & 117 117 tabim(nwl),tabimt(nwlt,4),tabre(nwl),tabret(nwlt,4),temref(4), & 118 118 wl(nwl),wlt(nwlt) … … 540 540 if(tk < temref(4)) tk=temref(4) 541 541 do 11 i=2,4 542 if(tk .ge.temref(i)) go to 12542 if(tk>=temref(i)) go to 12 543 543 11 continue 544 544 12 lt1=i 545 545 lt2=i-1 546 546 do 13 i=2,nwlt 547 if(alam .le.wlt(i)) go to 14547 if(alam<=wlt(i)) go to 14 548 548 13 continue 549 549 14 x1=log(wlt(i-1)) … … 586 586 Subroutine MieInt(Dx, SCm, Inp, Dqv, Dqxt, Dqsc, Dbsc, Dg, Xs1, Xs2, DPh, Error) 587 587 588 Integer * 2Imaxx588 Integer (kind=2) Imaxx 589 589 Parameter (Imaxx = 12000) 590 Real * 4RIMax ! largest real part of refractive index590 Real (kind=4) RIMax ! largest real part of refractive index 591 591 Parameter (RIMax = 2.5) 592 Real * 4IRIMax ! largest imaginary part of refractive index592 Real (kind=4) IRIMax ! largest imaginary part of refractive index 593 593 Parameter (IRIMax = -2) 594 Integer * 2Itermax594 Integer (kind=2) Itermax 595 595 Parameter (Itermax = 12000 * 2.5) 596 596 ! must be large enough to cope with the 597 597 ! largest possible nmx = x * abs(scm) + 15 598 598 ! or nmx = Dx + 4.05*Dx**(1./3.) + 2.0 599 Integer * 2Imaxnp599 Integer (kind=2) Imaxnp 600 600 Parameter (Imaxnp = 10000) ! Change this as required 601 601 ! INPUT 602 Real * 8Dx603 Complex * 16SCm604 Integer * 4Inp605 Real * 8Dqv(Inp)602 Real (kind=8) Dx 603 Complex (kind=8) SCm 604 Integer (kind=4) Inp 605 Real (kind=8) Dqv(Inp) 606 606 ! OUTPUT 607 Complex * 16Xs1(InP)608 Complex * 16Xs2(InP)609 Real * 8Dqxt610 Real * 8Dqsc611 Real * 8Dg612 Real * 8Dbsc613 Real * 8DPh(InP)614 Integer * 4Error607 Complex (kind=8) Xs1(InP) 608 Complex (kind=8) Xs2(InP) 609 Real (kind=8) Dqxt 610 Real (kind=8) Dqsc 611 Real (kind=8) Dg 612 Real (kind=8) Dbsc 613 Real (kind=8) DPh(InP) 614 Integer (kind=4) Error 615 615 ! LOCAL 616 Integer * 2I617 Integer * 2NStop618 Integer * 2NmX619 Integer * 4N ! N*N > 32767 ie N > 181620 Integer * 4Inp2621 Real * 8Chi,Chi0,Chi1622 Real * 8APsi,APsi0,APsi1623 Real * 8Pi0(Imaxnp)624 Real * 8Pi1(Imaxnp)625 Real * 8Taun(Imaxnp)626 Real * 8Psi,Psi0,Psi1627 Complex * 8Ir628 Complex * 16Cm629 Complex * 16A,ANM1,APB630 Complex * 16B,BNM1,AMB631 Complex * 16D(Itermax)632 Complex * 16Sp(Imaxnp)633 Complex * 16Sm(Imaxnp)634 Complex * 16Xi,Xi0,Xi1635 Complex * 16Y616 Integer (kind=2) I 617 Integer (kind=2) NStop 618 Integer (kind=2) NmX 619 Integer (kind=4) N ! N*N > 32767 ie N > 181 620 Integer (kind=4) Inp2 621 Real (kind=8) Chi,Chi0,Chi1 622 Real (kind=8) APsi,APsi0,APsi1 623 Real (kind=8) Pi0(Imaxnp) 624 Real (kind=8) Pi1(Imaxnp) 625 Real (kind=8) Taun(Imaxnp) 626 Real (kind=8) Psi,Psi0,Psi1 627 Complex (kind=4) Ir 628 Complex (kind=8) Cm 629 Complex (kind=8) A,ANM1,APB 630 Complex (kind=8) B,BNM1,AMB 631 Complex (kind=8) D(Itermax) 632 Complex (kind=8) Sp(Imaxnp) 633 Complex (kind=8) Sm(Imaxnp) 634 Complex (kind=8) Xi,Xi0,Xi1 635 Complex (kind=8) Y 636 636 ! ACCELERATOR VARIABLES 637 Integer * 2Tnp1638 Integer * 2Tnm1639 Real * 8Dn640 Real * 8Rnx641 Real * 8S(Imaxnp)642 Real * 8T(Imaxnp)643 Real * 8Turbo644 Real * 8A2645 Complex * 16A1637 Integer (kind=2) Tnp1 638 Integer (kind=2) Tnm1 639 Real (kind=8) Dn 640 Real (kind=8) Rnx 641 Real (kind=8) S(Imaxnp) 642 Real (kind=8) T(Imaxnp) 643 Real (kind=8) Turbo 644 Real (kind=8) A2 645 Complex (kind=8) A1 646 646 647 If ((Dx .Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then647 If ((Dx>Imaxx) .Or. (InP>ImaxNP)) Then 648 648 Error = 1 649 649 Return … … 652 652 Ir = 1 / Cm 653 653 Y = Dx * Cm 654 If (Dx .Lt.0.02) Then654 If (Dx<0.02) Then 655 655 NStop = 2 656 656 Else 657 If (Dx .Le.8.0) Then657 If (Dx<=8.0) Then 658 658 NStop = Dx + 4.00*Dx**(1./3.) + 2.0 659 659 Else 660 If (Dx .Lt.4200.0) Then660 If (Dx< 4200.0) Then 661 661 NStop = Dx + 4.05*Dx**(1./3.) + 2.0 662 662 Else … … 666 666 End If 667 667 NmX = Max(Real(NStop),Real(Abs(Y))) + 15. 668 If (Nmx .gt.Itermax) then668 If (Nmx > Itermax) then 669 669 Error = 1 670 670 Return … … 709 709 Dqxt = Tnp1 * Dble(A + B) + Dqxt 710 710 Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc 711 If (N .Gt.1) then711 If (N>1) then 712 712 Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN) 713 713 End If … … 717 717 AMB = A2 * (A - B) 718 718 Do I = 1,Inp2 719 If (I .GT.Inp) Then719 If (I>Inp) Then 720 720 S(I) = -Pi1(I) 721 721 Else … … 736 736 Xi1 = Dcmplx(APsi1,Chi1) 737 737 End Do 738 If (Dg .GT.0) Dg = 2 * Dg / Dqsc738 If (Dg >0) Dg = 2 * Dg / Dqsc 739 739 Dqsc = 2 * Dqsc / Dx**2 740 740 Dqxt = 2 * Dqxt / Dx**2 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cosp2/mrgrnk.F90
r3358 r5081 68 68 IRNGT (1) = 1 69 69 Return 70 Case Default71 Continue72 70 End Select 73 71 ! … … 268 266 IRNGT (1) = 1 269 267 Return 270 Case Default271 Continue272 268 End Select 273 269 ! … … 467 463 IRNGT (1) = 1 468 464 Return 469 Case Default470 Continue471 465 End Select 472 466 ! -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90
r3491 r5081 72 72 varOUT(1:dim1,1:dim2,1:dim3) = 0._wp 73 73 do j=1,dim2 74 where(flag(:,j,:) .eq.1)74 where(flag(:,j,:) == 1) 75 75 varOUT(:,j,:) = varIN2 76 76 endwhere 77 where(flag(:,j,:) .eq.2)77 where(flag(:,j,:) == 2) 78 78 varOUT(:,j,:) = varIN1 79 79 endwhere … … 96 96 97 97 varOUT(1:dim1,1:dim2,1:dim3) = 0._wp 98 where(flag(:,:,:) .eq.1)98 where(flag(:,:,:) == 1) 99 99 varOUT(:,:,:) = varIN2 100 100 endwhere 101 where(flag(:,:,:) .eq.2)101 where(flag(:,:,:) == 2) 102 102 varOUT(:,:,:) = varIN1 103 103 endwhere … … 295 295 296 296 ! Which LIDAR frequency are we using? 297 if (lidar_freq .eq.355) then297 if (lidar_freq == 355) then 298 298 Cmol = Cmol_355nm 299 299 rdiffm = rdiffm_355nm 300 300 endif 301 if (lidar_freq .eq.532) then301 if (lidar_freq == 532) then 302 302 Cmol = Cmol_532nm 303 303 rdiffm = rdiffm_532nm … … 336 336 337 337 ! LS and CONV Ice water coefficients 338 if (ice_type .eq.0) then338 if (ice_type == 0) then 339 339 polpart(INDX_LSICE,1:5) = polpartLSICE0 340 340 polpart(INDX_CVICE,1:5) = polpartCVICE0 341 341 endif 342 if (ice_type .eq.1) then342 if (ice_type == 1) then 343 343 polpart(INDX_LSICE,1:5) = polpartLSICE1 344 344 polpart(INDX_CVICE,1:5) = polpartCVICE1 … … 393 393 ! Polynomials kp_lidar derived from Mie theory 394 394 do i = 1, npart 395 where (rad_part(1:npoints,1:nlev,i) .gt.0.0)395 where (rad_part(1:npoints,1:nlev,i) > 0.0) 396 396 kp_part(1:npoints,1:nlev,i) = & 397 397 polpart(i,1)*(rad_part(1:npoints,1:nlev,i)*1e6)**4 & … … 426 426 ! Alpha of particles in each subcolumn: 427 427 do i = 1, npart 428 where (rad_part(1:npoints,1:nlev,i) .gt.0.0)428 where (rad_part(1:npoints,1:nlev,i) > 0.0) 429 429 alpha_part(1:npoints,1:nlev,i) = 3._wp/4._wp * Qscat & 430 430 * rhoair(1:npoints,1:nlev) * qpart(1:npoints,1:nlev,i) & -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/mrgrnk.F90
r3491 r5081 68 68 IRNGT (1) = 1 69 69 Return 70 Case Default71 Continue72 70 End Select 73 71 ! … … 268 266 IRNGT (1) = 1 269 267 Return 270 Case Default271 Continue272 268 End Select 273 269 ! … … 467 463 IRNGT (1) = 1 468 464 Return 469 Case Default470 Continue471 465 End Select 472 466 ! -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/optics_lib.F90
r3491 r5081 558 558 if(tk < temref(4)) tk=temref(4) 559 559 do i=2,4 560 if(tk .ge.temref(i)) go to 12560 if(tk>=temref(i)) go to 12 561 561 enddo 562 562 12 lt1 = i 563 563 lt2 = i-1 564 564 do i=2,nwlt 565 if(alam .le.wlt(i)) go to 14565 if(alam<=wlt(i)) go to 14 566 566 enddo 567 567 14 x1 = log(wlt(i-1)) … … 652 652 Complex(wp) :: A1 653 653 654 If ((Dx .Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then654 If ((Dx>Imaxx) .Or. (InP>ImaxNP)) Then 655 655 Error = 1 656 656 Return … … 659 659 Ir = 1 / Cm 660 660 Y = Dx * Cm 661 If (Dx .Lt.0.02) Then661 If (Dx<0.02) Then 662 662 NStop = 2 663 663 Else 664 If (Dx .Le.8.0) Then664 If (Dx<=8.0) Then 665 665 NStop = Dx + 4.00*Dx**(1./3.) + 2.0 666 666 Else 667 If (Dx .Lt.4200.0) Then667 If (Dx< 4200.0) Then 668 668 NStop = Dx + 4.05*Dx**(1./3.) + 2.0 669 669 Else … … 673 673 End If 674 674 NmX = Max(Real(NStop),Real(Abs(Y))) + 15. 675 If (Nmx .gt.Itermax) then675 If (Nmx > Itermax) then 676 676 Error = 1 677 677 Return … … 726 726 !ds Dqxt = Tnp1 * Dble(A + B) + Dqxt 727 727 Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc 728 If (N .Gt.1) then728 If (N>1) then 729 729 Dg = Dg + (dN*dN - 1) * (ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 *(ANM1*Conjg(BNM1)) / (dN*dN - dN) 730 730 !ds Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN) … … 735 735 AMB = A2 * (A - B) 736 736 Do I = 1,Inp2 737 If (I .GT.Inp) Then737 If (I>Inp) Then 738 738 S(I) = -Pi1(I) 739 739 Else … … 756 756 End Do 757 757 758 If (Dg .GT.0) Dg = 2 * Dg / Dqsc758 If (Dg >0) Dg = 2 * Dg / Dqsc 759 759 Dqsc = 2 * Dqsc / Dx**2 760 760 Dqxt = 2 * Dqxt / Dx**2 -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam.F90
r3491 r5081 179 179 180 180 ! Attenuation due to gaseous absorption between radar and volume 181 if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .and. pr .eq.1)) then181 if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .and. pr == 1)) then 182 182 if (d_gate==1) then 183 183 if (k>1) then … … 272 272 273 273 ! Which platforms to create diagnostics for? 274 if (platform .eq.'cloudsat') lcloudsat=.true.274 if (platform == 'cloudsat') lcloudsat=.true. 275 275 276 276 ! Create Cloudsat diagnostics. … … 289 289 enddo 290 290 enddo 291 where(cfad_ze .ne.R_UNDEF) cfad_ze = cfad_ze/Ncolumns291 where(cfad_ze /= R_UNDEF) cfad_ze = cfad_ze/Ncolumns 292 292 293 293 ! Compute cloudsat near-surface precipitation diagnostics … … 306 306 enddo 307 307 enddo 308 where(cfad_ze .ne.R_UNDEF) cfad_ze = cfad_ze/Ncolumns308 where(cfad_ze /= R_UNDEF) cfad_ze = cfad_ze/Ncolumns 309 309 endif 310 310 endif … … 402 402 do pr=1,Ncolumns 403 403 ! 1) Compute the PIA in all profiles containing hydrometeors 404 if ( (Ze_non_out(i,pr,cloudsat_preclvl) .gt.-100) .and. (Ze_out(i,pr,cloudsat_preclvl).gt.-100) ) then405 if ( (Ze_non_out(i,pr,cloudsat_preclvl) .lt.100) .and. (Ze_out(i,pr,cloudsat_preclvl).lt.100) ) then404 if ( (Ze_non_out(i,pr,cloudsat_preclvl)>-100) .and. (Ze_out(i,pr,cloudsat_preclvl)>-100) ) then 405 if ( (Ze_non_out(i,pr,cloudsat_preclvl)<100) .and. (Ze_out(i,pr,cloudsat_preclvl)<100) ) then 406 406 cloudsat_precip_pia(i,pr) = Ze_non_out(i,pr,cloudsat_preclvl) - Ze_out(i,pr,cloudsat_preclvl) 407 407 endif … … 412 412 ! 2a) Oceanic points. 413 413 ! ################################################################################ 414 if (land(i) .eq.0) then414 if (land(i) == 0) then 415 415 ! print*, 'aaa i, pr, fracPrecipIce(i,pr) : ', i, pr, fracPrecipIce(i,pr) !Artem 416 416 ! Snow 417 if(fracPrecipIce(i,pr) .gt.0.9) then418 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(2)) then417 if(fracPrecipIce(i,pr)>0.9) then 418 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(2)) then 419 419 cloudsat_pflag(i,pr) = pClass_Snow2 ! TSL: Snow certain 420 420 endif 421 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(4).and. &422 Ze_non_out(i,pr,cloudsat_preclvl) .le.Zenonbinval(2)) then421 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. & 422 Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(2)) then 423 423 cloudsat_pflag(i,pr) = pClass_Snow1 ! TSL: Snow possible 424 424 endif … … 426 426 427 427 ! Mixed 428 if(fracPrecipIce(i,pr) .gt.0.1.and.fracPrecipIce(i,pr).le.0.9) then429 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(2)) then428 if(fracPrecipIce(i,pr)>0.1.and.fracPrecipIce(i,pr)<=0.9) then 429 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(2)) then 430 430 cloudsat_pflag(i,pr) = pClass_Mixed2 ! TSL: Mixed certain 431 431 endif 432 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(4).and. &433 Ze_non_out(i,pr,cloudsat_preclvl) .le.Zenonbinval(2)) then432 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. & 433 Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(2)) then 434 434 cloudsat_pflag(i,pr) = pClass_Mixed1 ! TSL: Mixed possible 435 435 endif … … 437 437 438 438 ! Rain 439 if(fracPrecipIce(i,pr) .le.0.1) then440 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(1)) then439 if(fracPrecipIce(i,pr)<=0.1) then 440 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(1)) then 441 441 cloudsat_pflag(i,pr) = pClass_Rain3 ! TSL: Rain certain 442 442 endif 443 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(3).and. &444 Ze_non_out(i,pr,cloudsat_preclvl) .le.Zenonbinval(1)) then443 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(3).and. & 444 Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(1)) then 445 445 cloudsat_pflag(i,pr) = pClass_Rain2 ! TSL: Rain probable 446 446 endif 447 if(Ze_non_out(i,pr,cloudsat_preclvl) .gt.Zenonbinval(4).and. &448 Ze_non_out(i,pr,cloudsat_preclvl) .le.Zenonbinval(3)) then447 if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. & 448 Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(3)) then 449 449 cloudsat_pflag(i,pr) = pClass_Rain1 ! TSL: Rain possible 450 450 endif 451 if(cloudsat_precip_pia(i,pr) .gt.40) then451 if(cloudsat_precip_pia(i,pr)>40) then 452 452 cloudsat_pflag(i,pr) = pClass_Rain4 ! TSL: Heavy Rain 453 453 endif … … 455 455 456 456 ! No precipitation 457 if(Ze_non_out(i,pr,cloudsat_preclvl) .le.-15) then457 if(Ze_non_out(i,pr,cloudsat_preclvl)<=-15) then 458 458 cloudsat_pflag(i,pr) = pClass_noPrecip ! TSL: Not Raining 459 459 endif … … 463 463 ! 2b) Land points. 464 464 ! ################################################################################ 465 if (land(i) .eq.1) then465 if (land(i) == 1) then 466 466 ! Find Zmax, the maximum reflectivity value in the attenuated profile (Ze_out); 467 467 Zmax=maxval(Ze_out(i,pr,:)) 468 468 469 469 ! Snow (T<273) 470 if(t2m(i) .lt.273._wp) then471 if(Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(5)) then470 if(t2m(i) < 273._wp) then 471 if(Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(5)) then 472 472 cloudsat_pflag(i,pr) = pClass_Snow2 ! JEK: Snow certain 473 473 endif 474 if(Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(6) .and. &475 Ze_out(i,pr,cloudsat_preclvl) .le.Zbinvallnd(5)) then474 if(Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6) .and. & 475 Ze_out(i,pr,cloudsat_preclvl)<=Zbinvallnd(5)) then 476 476 cloudsat_pflag(i,pr) = pClass_Snow1 ! JEK: Snow possible 477 477 endif … … 479 479 480 480 ! Mized phase (273<T<275) 481 if(t2m(i) .ge. 273._wp .and. t2m(i) .le.275._wp) then482 if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. &483 (Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(4))) then481 if(t2m(i) >= 273._wp .and. t2m(i) <= 275._wp) then 482 if ((Zmax > Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr)>30) .or. & 483 (Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(4))) then 484 484 cloudsat_pflag(i,pr) = pClass_Mixed2 ! JEK: Mixed certain 485 485 endif 486 if ((Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(6) .and. &487 Ze_out(i,pr,cloudsat_preclvl) .le.Zbinvallnd(4)) .and. &488 (Zmax .gt.Zbinvallnd(5)) ) then486 if ((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6) .and. & 487 Ze_out(i,pr,cloudsat_preclvl) <= Zbinvallnd(4)) .and. & 488 (Zmax > Zbinvallnd(5)) ) then 489 489 cloudsat_pflag(i,pr) = pClass_Mixed1 ! JEK: Mixed possible 490 490 endif … … 492 492 493 493 ! Rain (T>275) 494 if(t2m(i) .gt.275) then495 if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. &496 (Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(2))) then494 if(t2m(i) > 275) then 495 if ((Zmax > Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr)>30) .or. & 496 (Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(2))) then 497 497 cloudsat_pflag(i,pr) = pClass_Rain3 ! JEK: Rain certain 498 498 endif 499 if((Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(6)) .and. &500 (Zmax .gt.Zbinvallnd(3))) then499 if((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6)) .and. & 500 (Zmax > Zbinvallnd(3))) then 501 501 cloudsat_pflag(i,pr) = pClass_Rain2 ! JEK: Rain probable 502 502 endif 503 if((Ze_out(i,pr,cloudsat_preclvl) .gt.Zbinvallnd(6)) .and. &504 (Zmax .lt.Zbinvallnd(3))) then503 if((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6)) .and. & 504 (Zmax<Zbinvallnd(3))) then 505 505 cloudsat_pflag(i,pr) = pClass_Rain1 ! JEK: Rain possible 506 506 endif 507 if(cloudsat_precip_pia(i,pr) .gt.40) then507 if(cloudsat_precip_pia(i,pr)>40) then 508 508 cloudsat_pflag(i,pr) = pClass_Rain4 ! JEK: Heavy Rain 509 509 endif … … 511 511 512 512 ! No precipitation 513 if(Ze_out(i,pr,cloudsat_preclvl) .le.-15) then513 if(Ze_out(i,pr,cloudsat_preclvl)<=-15) then 514 514 cloudsat_pflag(i,pr) = pClass_noPrecip ! JEK: Not Precipitating 515 515 endif … … 526 526 ! Gridmean precipitation fraction for each precipitation type 527 527 do k=1,nCloudsatPrecipClass 528 if (any(cloudsat_pflag(i,:) .eq.k-1)) then529 cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) .eq.k-1)528 if (any(cloudsat_pflag(i,:) == k-1)) then 529 cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) == k-1) 530 530 endif 531 531 enddo -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam_optics.F90
r3491 r5081 172 172 173 173 ! Compute effective radius from number concentration and distribution parameters 174 if (Re_internal .eq.0) then174 if (Re_internal == 0) then 175 175 call calc_Re(hm_matrix(pr,k,tp),Np_matrix(pr,k,tp),rho_a, & 176 176 sd%dtype(tp),sd%apm(tp),sd%bpm(tp),sd%rho(tp),sd%p1(tp),sd%p2(tp),sd%p3(tp),Re) … … 187 187 ! Index into particle size dimension of scaling tables 188 188 iRe_type=1 189 if(Re .gt.0) then189 if(Re>0) then 190 190 ! Determine index in to scale LUT 191 191 ! Distance between Re points (defined by "base" and "step") for … … 197 197 base = rcfg%base_list(n+1) 198 198 iRe_type=Re/step 199 if (iRe_type .lt.1) iRe_type=1199 if (iRe_type<1) iRe_type=1 200 200 Re=step*(iRe_type+0.5_wp) ! set value of Re to closest value allowed in LUT. 201 201 iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step) 202 202 203 203 ! Make sure iRe_type is within bounds 204 if (iRe_type .ge.nRe_types) then204 if (iRe_type>=nRe_types) then 205 205 !write(*,*) 'Warning: size of Re exceed value permitted ', & 206 206 ! 'in Look-Up Table (LUT). Will calculate. ' … … 405 405 ! Exponential is same as modified gamma with vu =1 406 406 ! if Np is specified then we will just treat as modified gamma 407 if(dtype .eq. 2 .and. Np .gt.0) then407 if(dtype == 2 .and. Np > 0) then 408 408 local_dtype = 1 409 409 local_p3 = 1 … … 441 441 endif 442 442 443 if( Np .eq.0 .and. p2+1 > 1E-8) then ! use default value for MEAN diameter as first default443 if( Np==0 .and. p2+1 > 1E-8) then ! use default value for MEAN diameter as first default 444 444 dm = p2 ! by definition, should have units of microns 445 445 D0 = gamma(vu)/gamma(vu+1)*dm 446 446 else ! use value of Np 447 if(Np .eq.0) then447 if(Np==0) then 448 448 if( abs(p1+1) > 1E-8 ) then ! use default number concentration 449 449 local_Np = p1 ! total number concentration / pa --- units kg^-1 … … 525 525 526 526 ! get rg ... 527 if( Np .eq.0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg527 if( Np==0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg 528 528 rg = p2 529 529 else … … 826 826 log_sigma_g = p3 827 827 tmp2 = (bpm*log_sigma_g)*(bpm*log_sigma_g) 828 if(Re .le.0) then828 if(Re<=0) then 829 829 rg = p2 830 830 else -
LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90
r3491 r5081 75 75 76 76 ! Test for valid input overlap assumption 77 if (overlap .ne. 1 .and. overlap .ne. 2 .and. overlap .ne.3) then77 if (overlap /= 1 .and. overlap /= 2 .and. overlap /= 3) then 78 78 overlap=default_overlap 79 79 call errorMessage('ERROR(scops): Invalid overlap assumption provided. Using default overlap assumption (max/ran)') … … 92 92 tca(1:npoints,1:nlev) = cc(1:npoints,1:nlev) 93 93 94 if (ncolprint .ne.0) then94 if (ncolprint/=0) then 95 95 write (6,'(a)') 'frac_out_pp_rev:' 96 96 do j=1,npoints,1000 … … 102 102 write (6,'(I3)') ncol 103 103 endif 104 if (ncolprint .ne.0) then104 if (ncolprint/=0) then 105 105 write (6,'(a)') 'last_frac_pp:' 106 106 do j=1,npoints,1000 … … 122 122 123 123 ! Initialise threshold 124 IF (ilev .eq.1) then124 IF (ilev==1) then 125 125 ! If max overlap 126 IF (overlap .eq.1) then126 IF (overlap==1) then 127 127 ! Select pixels spread evenly across the gridbox 128 128 threshold(1:npoints,1:ncol)=boxpos(1:npoints,1:ncol) … … 137 137 enddo 138 138 ENDIF 139 IF (ncolprint .ne.0) then139 IF (ncolprint/=0) then 140 140 write (6,'(a)') 'threshold_nsf2:' 141 141 do j=1,npoints,1000 … … 147 147 ENDIF 148 148 149 IF (ncolprint .ne.0) then149 IF (ncolprint/=0) then 150 150 write (6,'(a)') 'ilev:' 151 151 write (6,'(I2)') ilev … … 157 157 !maxocc(1:npoints,ibox) = merge(1,0, conv(1:npoints,ilev) .gt. boxpos(1:npoints,ibox)) 158 158 do j=1,npoints 159 if (boxpos(j,ibox) .le.conv(j,ilev)) then159 if (boxpos(j,ibox)<=conv(j,ilev)) then 160 160 maxocc(j,ibox) = 1 161 161 else … … 165 165 166 166 ! Max overlap 167 if (overlap .eq.1) then167 if (overlap==1) then 168 168 threshold_min(1:npoints,ibox) = conv(1:npoints,ilev) 169 169 maxosc(1:npoints,ibox) = 1 … … 171 171 172 172 ! Random overlap 173 if (overlap .eq.2) then173 if (overlap==2) then 174 174 threshold_min(1:npoints,ibox) = conv(1:npoints,ilev) 175 175 maxosc(1:npoints,ibox) = 0 176 176 endif 177 177 ! Max/Random overlap 178 if (overlap .eq.3) then178 if (overlap==3) then 179 179 ! DS2014 START: The bounds on tca are not valid when ilev=1. 180 180 !threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev))) … … 182 182 ! min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .and. & 183 183 ! (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev))) 184 if (ilev .ne.1) then184 if (ilev /= 1) then 185 185 threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev))) 186 maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt.&186 maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) < & 187 187 min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .and. & 188 (threshold(1:npoints,ibox) .gt.conv(1:npoints,ilev)))188 (threshold(1:npoints,ibox)>conv(1:npoints,ilev))) 189 189 else 190 190 threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(0._wp,tca(1:npoints,ilev))) 191 maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt.&191 maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) < & 192 192 min(0._wp,tca(1:npoints,ilev)) .and. & 193 (threshold(1:npoints,ibox) .gt.conv(1:npoints,ilev)))193 (threshold(1:npoints,ibox)>conv(1:npoints,ilev))) 194 194 endif 195 195 endif … … 205 205 206 206 ! Fill frac_out with 1's where tca is greater than the threshold 207 frac_out(1:npoints,ibox,ilev) = merge(1,0,tca(1:npoints,ilev) .gt.threshold(1:npoints,ibox))207 frac_out(1:npoints,ibox,ilev) = merge(1,0,tca(1:npoints,ilev)>threshold(1:npoints,ibox)) 208 208 209 209 ! Code to partition boxes into startiform and convective parts goes here 210 where(threshold(1:npoints,ibox) .le.conv(1:npoints,ilev) .and. conv(1:npoints,ilev).gt.0.) frac_out(1:npoints,ibox,ilev)=2210 where(threshold(1:npoints,ibox)<=conv(1:npoints,ilev) .and. conv(1:npoints,ilev)>0.) frac_out(1:npoints,ibox,ilev)=2 211 211 ENDDO ! ibox 212 212 213 213 214 214 ! Set last_frac to tca at this level, so as to be tca from last level next time round 215 if (ncolprint .ne.0) then215 if (ncolprint/=0) then 216 216 do j=1,npoints ,1000 217 217 write(6,'(a10)') 'j=' -
LMDZ6/branches/Amaury_dev/libf/phylmd/cv30param.h
r1992 r5081 20 20 real betad 21 21 22 COMMON /cv30param/ noff, minorig, nl, nlp, nlm & 23 , sigd, spfac & 22 COMMON /cv30param/ sigd, spfac & 24 23 !IM cf. FH : pour compatibilite avec conema3 TEMPORAIRE : ,pbcrit, ptcrit, epmax 25 24 ,pbcrit, ptcrit & 26 25 ,omtrain & 27 26 ,dtovsh, dpbase, dttrig & 28 ,dtcrit, tau, beta, alpha, delta, betad 27 ,dtcrit, tau, beta, alpha, delta, betad & 28 ,noff, minorig, nl, nlp, nlm 29 29 30 30 !$OMP THREADPRIVATE(/cv30param/) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cv3param.h
r3571 r5081 38 38 ,delta, betad & 39 39 ,ejectliq, ejectice & 40 ,flag_wb & 40 41 ,flag_epKEorig & 41 , flag_wb,cv_flag_feed &42 ,cv_flag_feed & 42 43 ,noff, minorig, nl, nlp, nlm 43 44 !$OMP THREADPRIVATE(/cv3param/) -
LMDZ6/branches/Amaury_dev/libf/phylmd/cva_driver.F90
r4613 r5081 858 858 ! is assumed useless. 859 859 ! 860 compress = ncum .lt.len*comp_threshold860 compress = ncum < len*comp_threshold 861 861 ! 862 862 IF (.not. compress) THEN … … 1100 1100 Print *, 'cva_driver after cv3_unsat:mp , water, ice, evap, fondue ' 1101 1101 DO k = 1,nd 1102 write (6, '(i4,5(1x,e13.6))') ,&1102 write (6, '(i4,5(1x,e13.6))') & 1103 1103 k, mp(igout,k), water(igout,k), ice(igout,k), & 1104 1104 evap(igout,k), fondue(igout,k) … … 1106 1106 Print *, 'cva_driver after cv3_unsat: wdtrainA, wdtrainS, wdtrainM ' !!jygprl 1107 1107 DO k = 1,nd 1108 write (6, '(i4,3(1x,e13.6))') ,&1108 write (6, '(i4,3(1x,e13.6))') & 1109 1109 k, wdtrainA(igout,k), wdtrainS(igout,k), wdtrainM(igout,k) !!jygprl 1110 1110 ENDDO -
LMDZ6/branches/Amaury_dev/libf/phylmd/cvparam.h
r1992 r5081 21 21 real delta 22 22 23 COMMON /cvparam/ noff, minorig, nl, nlp, nlm & 24 ,elcrit, tlcrit & 23 COMMON /cvparam/ elcrit, tlcrit & 25 24 ,entp, sigs, sigd & 26 25 ,omtrain, omtsnow, coeffr, coeffs & 27 ,dtmax, cu, betad, alpha, damp, delta 26 ,dtmax, cu, betad, alpha, damp, delta & 27 ,noff, minorig, nl, nlp, nlm 28 28 29 29 !$OMP THREADPRIVATE(/cvparam/) -
LMDZ6/branches/Amaury_dev/libf/phylmd/flux_arp.h
r3780 r5081 15 15 real :: tg 16 16 17 common /flux_arp/fsens,flat, ust,tg,ok_flux_surf,ok_prescr_ust,ok_prescr_beta,betaevap,ok_forc_tsurf17 common /flux_arp/fsens,flat,betaevap,ust,tg,ok_flux_surf,ok_prescr_ust,ok_prescr_beta,ok_forc_tsurf 18 18 19 19 !$OMP THREADPRIVATE(/flux_arp/) -
LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_ratqs_ini.F90
r4812 r5081 3 3 IMPLICIT NONE 4 4 5 save 5 INTEGER :: lunout 6 6 7 integer :: lunout 7 INTEGER, PROTECTED :: nbsrf,is_lic,is_ter 8 REAL, PROTECTED :: RG,RV,RD,RCPD,RLSTT,RLVTT,RTT 9 REAL, PROTECTED :: a_ratqs_cv 10 REAL, PROTECTED :: tau_var 11 REAL, PROTECTED :: fac_tau 12 REAL, PROTECTED :: tau_cumul 13 REAL, PROTECTED :: a_ratqs_wake 14 INTEGER, PROTECTED :: dqimpl 8 15 9 INTEGER, SAVE, PROTECTED :: nbsrf,is_lic,is_ter 10 REAL, SAVE, PROTECTED :: RG,RV,RD,RCPD,RLSTT,RLVTT,RTT 11 REAL, SAVE, PROTECTED :: a_ratqs_cv 12 REAL, SAVE, PROTECTED :: tau_var 13 REAL, SAVE, PROTECTED :: fac_tau 14 REAL, SAVE, PROTECTED :: tau_cumul 15 REAL, SAVE, PROTECTED :: a_ratqs_wake 16 INTEGER, SAVE, PROTECTED :: dqimpl 17 18 real, allocatable, SAVE :: povariance(:,:) 16 REAL, ALLOCATABLE :: povariance(:,:) 19 17 !$OMP THREADPRIVATE(povariance) 20 real, allocatable, SAVE :: var_conv(:,:)18 REAL, ALLOCATABLE :: var_conv(:,:) 21 19 !$OMP THREADPRIVATE(var_conv) 22 20 … … 78 76 !-------------------------------------------------------- 79 77 80 if (klon .eq.1) then78 if (klon==1) then 81 79 do k=1,klev 82 80 do i=1,klon -
LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_ini.F90
r4863 r5081 2 2 3 3 IMPLICIT NONE 4 5 save6 7 4 8 5 integer, protected :: dvdq=1,dqimpl=-1,prt_level=0,lunout -
LMDZ6/branches/Amaury_dev/libf/phylmd/rrtm/rrtm_taumol1.F90
r2003 r5081 190 190 IF (LHOOK) CALL DR_HOOK('RRTM_TAUMOL1',0,ZHOOK_HANDLE) 191 191 !--ajout OB 192 IF (K_LAYTROP .GT.100) THEN192 IF (K_LAYTROP>100) THEN 193 193 PRINT *,'ATTENTION KLAY_TROP > 100 PROBLEME ARRAY DANS RRTM ON ARRETE' 194 194 STOP -
LMDZ6/branches/Amaury_dev/libf/phylmd/tlift.F90
r2197 r5081 213 213 tvp(i) = tpk(i)*(1.+qsat_new/eps-rr(nk)) 214 214 ! jyg2 215 ELSE216 CONTINUE217 215 END IF 218 216 -
LMDZ6/branches/Amaury_dev/libf/phylmd/yamada_ini_mod.F90
r4822 r5081 5 5 6 6 implicit none 7 8 save9 7 10 8 LOGICAL :: new_yamada4
Note: See TracChangeset
for help on using the changeset viewer.