diff --git a/src/arome/ext/aroini_micro.F90 b/src/arome/ext/aroini_micro.F90 index 4881f58e8cf6e2881811f1d64ce3cc3c2764572d..2c58dac7878f894c1c73ff099481792c66ab91e9 100644 --- a/src/arome/ext/aroini_micro.F90 +++ b/src/arome/ext/aroini_micro.F90 @@ -150,6 +150,7 @@ LSEDIC=LDSEDIC ! et MODD_RAIN_ICE_PARAM LSNOW_T=.FALSE. LRED=CMICRO=='ICE3' .OR. CMICRO=='ICE4' +LPACK_INTERP=.TRUE. CALL INI_RAIN_ICE (KULOUT, PTSTEP, 20.,KSPLITR,CMICRO) CALL INI_TIWMX diff --git a/src/common/micro/interp_micro.func.h b/src/common/micro/interp_micro.func.h new file mode 100644 index 0000000000000000000000000000000000000000..ec9a912eafdbc59d879bea36596c0f73461e4ecb --- /dev/null +++ b/src/common/micro/interp_micro.func.h @@ -0,0 +1,269 @@ +!These routines are intented to be included in the contains part of other subroutines. +!To allow the transformation for GPU, no local array must be declared. +!If a temporary local array is needed, it must be added as a buffer in the interface (IBUF?, ZBUF?) + +SUBROUTINE INTERP_MICRO_1D(KPROMA, KSIZE, PIN, KNUM, P1, P2, & + LDPACK, LDMASK, KBUF1, KBUF2, PBUF1, PBUF2, & + KLEN, & + PLT1, POUT1, PLT2, POUT2, PLT3, POUT3) + +IMPLICIT NONE + +INTEGER, INTENT(IN) :: KPROMA !Array size +INTEGER, INTENT(IN) :: KSIZE !Last usefull array index +REAL, DIMENSION(KPROMA), INTENT(IN) :: PIN !Input array +INTEGER, INTENT(IN) :: KNUM !Number of points in the look-up table +REAL, INTENT(IN) :: P1 !Scaling factor +REAL, INTENT(IN) :: P2 !Scaling factor +LOGICAL, INTENT(IN) :: LDPACK !.TRUE. to perform packing +LOGICAL, DIMENSION(KPROMA), INTENT(IN) :: LDMASK !Computation mask +INTEGER, DIMENSION(KPROMA), INTENT(OUT) :: KBUF1, KBUF2 !Buffer arrays +REAL, DIMENSION(KPROMA), INTENT(OUT) :: PBUF1, PBUF2 !Buffer arrays +INTEGER, INTENT(OUT) :: KLEN !Number of active points +REAL, DIMENSION(KNUM), INTENT(IN) :: PLT1 !Look-up table +REAL, DIMENSION(KPROMA), INTENT(OUT) :: POUT1 !Interpolated values +REAL, DIMENSION(KNUM), INTENT(IN) , OPTIONAL :: PLT2 +REAL, DIMENSION(KPROMA), INTENT(OUT), OPTIONAL :: POUT2 +REAL, DIMENSION(KNUM), INTENT(IN) , OPTIONAL :: PLT3 +REAL, DIMENSION(KPROMA), INTENT(OUT), OPTIONAL :: POUT3 + +INTEGER :: JL +INTEGER :: IINDEX +REAL :: ZINDEX + +IF (LDPACK) THEN + + !Pack input array + KLEN=0 + DO JL=1, KSIZE + IF (LDMASK(JL)) THEN + KLEN=KLEN+1 + PBUF1(KLEN)=PIN(JL) + KBUF1(KLEN)=JL + ENDIF + ENDDO + + IF (KLEN>0) THEN + !Index computation + !$mnh_expand_array(JL=1:KLEN) + PBUF1(1:KLEN) = MAX(1.00001, MIN(REAL(KNUM)-0.00001, P1 * LOG(PBUF1(1:KLEN)) + P2)) + KBUF2(1:KLEN) = INT(PBUF1(1:KLEN)) + PBUF1(1:KLEN) = PBUF1(1:KLEN) - REAL(KBUF2(1:KLEN)) + !$mnh_end_expand_array(JL=1:KLEN) + + !Interpolation and unpack + !$mnh_expand_array(JL=1:KLEN) + PBUF2(1:KLEN) = PLT1(KBUF2(1:KLEN)+1) * PBUF1(1:KLEN) & + &-PLT1(KBUF2(1:KLEN) ) * (PBUF1(1:KLEN) - 1.0) + !$mnh_end_expand_array(JL=1:KLEN) + POUT1(:)=0. + DO JL=1, KLEN + POUT1(KBUF1(JL))=PBUF2(JL) + ENDDO + + !Interpolation and unpack 2 + IF(PRESENT(PLT2)) THEN + !$mnh_expand_array(JL=1:KLEN) + PBUF2(1:KLEN) = PLT2(KBUF2(1:KLEN)+1) * PBUF1(1:KLEN) & + &-PLT2(KBUF2(1:KLEN) ) * (PBUF1(1:KLEN) - 1.0) + !$mnh_end_expand_array(JL=1:KLEN) + POUT2(:)=0. + DO JL=1, KLEN + POUT2(KBUF1(JL))=PBUF2(JL) + ENDDO + ENDIF + + !Interpolation and unpack 3 + IF(PRESENT(PLT3)) THEN + !$mnh_expand_array(JL=1:KLEN) + PBUF2(1:KLEN) = PLT3(KBUF2(1:KLEN)+1) * PBUF1(1:KLEN) & + &-PLT3(KBUF2(1:KLEN) ) * (PBUF1(1:KLEN) - 1.0) + !$mnh_end_expand_array(JL=1:KLEN) + POUT3(:)=0. + DO JL=1, KLEN + POUT3(KBUF1(JL))=PBUF2(JL) + ENDDO + ENDIF + + ENDIF + +ELSE + + KLEN=0 + DO JL=1, KSIZE + IF (LDMASK(JL)) THEN + KLEN=KLEN+1 + + !Index computation + ZINDEX = MAX(1.00001, MIN(REAL(KNUM)-0.00001, P1 * LOG(PIN(JL)) + P2)) + IINDEX = INT(ZINDEX) + ZINDEX = ZINDEX - REAL(IINDEX) + + !Interpolations + POUT1(JL) = PLT1(IINDEX+1) * ZINDEX & + &-PLT1(IINDEX ) * (ZINDEX - 1.0) + + IF(PRESENT(PLT2)) THEN + POUT2(JL) = PLT2(IINDEX+1) * ZINDEX & + &-PLT2(IINDEX ) * (ZINDEX - 1.0) + ENDIF + + IF(PRESENT(PLT3)) THEN + POUT3(JL) = PLT3(IINDEX+1) * ZINDEX & + &-PLT3(IINDEX ) * (ZINDEX - 1.0) + ENDIF + + ELSE + POUT1(JL) = 0. + IF(PRESENT(PLT2)) POUT2(JL) = 0. + IF(PRESENT(PLT3)) POUT3(JL) = 0. + ENDIF + ENDDO + +ENDIF +END SUBROUTINE INTERP_MICRO_1D + +SUBROUTINE INTERP_MICRO_2D(KPROMA, KSIZE, PIN1, PIN2, KNUM1, KNUM2, P11, P12, P21, P22,& + LDPACK, LDMASK, KBUF1, KBUF2, KBUF3, PBUF1, PBUF2, PBUF3, & + KLEN, & + PLT1, POUT1, PLT2, POUT2, PLT3, POUT3) + +IMPLICIT NONE + +INTEGER, INTENT(IN) :: KPROMA !Array size +INTEGER, INTENT(IN) :: KSIZE !Last usefull array index +REAL, DIMENSION(KPROMA), INTENT(IN) :: PIN1 !Input array +REAL, DIMENSION(KPROMA), INTENT(IN) :: PIN2 !Input array +INTEGER, INTENT(IN) :: KNUM1 !First dimension of the look-up table +INTEGER, INTENT(IN) :: KNUM2 !Second dimension of the look-up table +REAL, INTENT(IN) :: P11 !Scaling factor +REAL, INTENT(IN) :: P12 !Scaling factor +REAL, INTENT(IN) :: P21 !Scaling factor +REAL, INTENT(IN) :: P22 !Scaling factor +LOGICAL, INTENT(IN) :: LDPACK !.TRUE. to perform packing +LOGICAL, DIMENSION(KPROMA), INTENT(IN) :: LDMASK !Computation mask +INTEGER, DIMENSION(KPROMA), INTENT(OUT) :: KBUF1, KBUF2, KBUF3 !Buffer arrays +REAL, DIMENSION(KPROMA), INTENT(OUT) :: PBUF1, PBUF2, PBUF3 !Buffer arrays +INTEGER, INTENT(OUT) :: KLEN !Number of active points +REAL, DIMENSION(KNUM1, KNUM2), INTENT(IN) :: PLT1 !Look-up table +REAL, DIMENSION(KPROMA), INTENT(OUT) :: POUT1 !Interpolated values from the first look-up table +REAL, DIMENSION(KNUM1, KNUM2), INTENT(IN) , OPTIONAL :: PLT2 !Other look-up table +REAL, DIMENSION(KPROMA), INTENT(OUT), OPTIONAL :: POUT2 !Interpolated values from the second look-up table +REAL, DIMENSION(KNUM2, KNUM1), INTENT(IN) , OPTIONAL :: PLT3 !Another look-up table **CAUTION, TABLE IS REVERSED** +REAL, DIMENSION(KPROMA), INTENT(OUT), OPTIONAL :: POUT3 !Interpolated values from the third look-up table + +INTEGER :: JL +INTEGER :: IINDEX1, IINDEX2 +REAL :: ZINDEX1, ZINDEX2 + +IF (LDPACK) THEN + + !Pack input array + KLEN=0 + DO JL=1, KSIZE + IF (LDMASK(JL)) THEN + KLEN=KLEN+1 + PBUF1(KLEN)=PIN1(JL) + PBUF2(KLEN)=PIN2(JL) + KBUF3(KLEN)=JL + ENDIF + ENDDO + + IF (KLEN>0) THEN + !Index computation + !$mnh_expand_array(JL=1:KLEN) + PBUF1(1:KLEN) = MAX(1.00001, MIN(REAL(KNUM1)-0.00001, P11 * LOG(PBUF1(1:KLEN)) + P12)) + KBUF1(1:KLEN) = INT(PBUF1(1:KLEN)) + PBUF1(1:KLEN) = PBUF1(1:KLEN) - REAL(KBUF1(1:KLEN)) + + PBUF2(1:KLEN) = MAX(1.00001, MIN(REAL(KNUM2)-0.00001, P21 * LOG(PBUF2(1:KLEN)) + P22)) + KBUF2(1:KLEN) = INT(PBUF2(1:KLEN)) + PBUF2(1:KLEN) = PBUF2(1:KLEN) - REAL(KBUF2(1:KLEN)) + !$mnh_end_expand_array(JL=1:KLEN) + + !Interpolation and unpack 1 + DO JL=1, KLEN + PBUF3(JL) = ( PLT1(KBUF1(JL)+1,KBUF2(JL)+1)* PBUF2(JL) & + -PLT1(KBUF1(JL)+1,KBUF2(JL) )*(PBUF2(JL) - 1.0)) * PBUF1(JL) & + -( PLT1(KBUF1(JL) ,KBUF2(JL)+1)* PBUF2(JL) & + -PLT1(KBUF1(JL) ,KBUF2(JL) )*(PBUF2(JL) - 1.0)) * (PBUF1(JL) - 1.0) + ENDDO + POUT1(:)=0. + DO JL=1, KLEN + POUT1(KBUF3(JL))=PBUF3(JL) + ENDDO + + !Interpolation and unpack 2 + IF(PRESENT(PLT2)) THEN + DO JL=1, KLEN + PBUF3(JL) = ( PLT2(KBUF1(JL)+1,KBUF2(JL)+1)* PBUF2(JL) & + -PLT2(KBUF1(JL)+1,KBUF2(JL) )*(PBUF2(JL) - 1.0)) * PBUF1(JL) & + -( PLT2(KBUF1(JL) ,KBUF2(JL)+1)* PBUF2(JL) & + -PLT2(KBUF1(JL) ,KBUF2(JL) )*(PBUF2(JL) - 1.0)) * (PBUF1(JL) - 1.0) + ENDDO + POUT2(:)=0. + DO JL=1, KLEN + POUT2(KBUF3(JL))=PBUF3(JL) + ENDDO + ENDIF + + !Interpolation and unpack 3 + IF(PRESENT(PLT3)) THEN + DO JL=1, KLEN + PBUF3(JL) = ( PLT3(KBUF2(JL)+1,KBUF1(JL)+1)* PBUF1(JL) & + -PLT3(KBUF2(JL)+1,KBUF1(JL) )*(PBUF1(JL) - 1.0)) * PBUF2(JL) & + -( PLT3(KBUF2(JL) ,KBUF1(JL)+1)* PBUF1(JL) & + -PLT3(KBUF2(JL) ,KBUF1(JL) )*(PBUF1(JL) - 1.0)) * (PBUF2(JL) - 1.0) + ENDDO + POUT3(:)=0. + DO JL=1, KLEN + POUT3(KBUF3(JL))=PBUF3(JL) + ENDDO + ENDIF + ENDIF + +ELSE + + KLEN=0 + DO JL=1, KSIZE + IF (LDMASK(JL)) THEN + KLEN=KLEN+1 + + !Indexes computation + ZINDEX1 = MAX(1.00001, MIN(REAL(KNUM1)-0.00001, P11 * LOG(PIN1(JL)) + P12)) + IINDEX1 = INT(ZINDEX1) + ZINDEX1 = ZINDEX1 - REAL(IINDEX1) + + ZINDEX2 = MAX(1.00001, MIN(REAL(KNUM1)-0.00001, P21 * LOG(PIN2(JL)) + P22)) + IINDEX2 = INT(ZINDEX2) + ZINDEX2 = ZINDEX2 - REAL(IINDEX2) + + !Interpolations + POUT1(JL) = ( PLT1(IINDEX1+1,IINDEX2+1)* ZINDEX2 & + -PLT1(IINDEX1+1,IINDEX2 )*(ZINDEX2 - 1.0)) * ZINDEX1 & + -( PLT1(IINDEX1 ,IINDEX2+1)* ZINDEX2 & + -PLT1(IINDEX1 ,IINDEX2 )*(ZINDEX2 - 1.0)) * (ZINDEX1 - 1.0) + + IF(PRESENT(PLT2)) THEN + POUT2(JL) = ( PLT2(IINDEX1+1,IINDEX2+1)* ZINDEX2 & + -PLT2(IINDEX1+1,IINDEX2 )*(ZINDEX2 - 1.0)) * ZINDEX1 & + -( PLT2(IINDEX1 ,IINDEX2+1)* ZINDEX2 & + -PLT2(IINDEX1 ,IINDEX2 )*(ZINDEX2 - 1.0)) * (ZINDEX1 - 1.0) + ENDIF + + IF(PRESENT(PLT3)) THEN + POUT3(JL) = ( PLT3(IINDEX2+1,IINDEX1+1)* ZINDEX1 & + -PLT3(IINDEX2+1,IINDEX1 )*(ZINDEX1 - 1.0)) * ZINDEX2 & + -( PLT3(IINDEX2 ,IINDEX1+1)* ZINDEX1 & + -PLT3(IINDEX2 ,IINDEX1 )*(ZINDEX1 - 1.0)) * (ZINDEX2 - 1.0) + ENDIF + + ELSE + POUT1(JL)=0. + IF(PRESENT(PLT2)) POUT2(JL)=0. + IF(PRESENT(PLT3)) POUT3(JL)=0. + ENDIF + ENDDO + +ENDIF +END SUBROUTINE INTERP_MICRO_2D diff --git a/src/common/micro/modd_param_ice.F90 b/src/common/micro/modd_param_ice.F90 index 0ce2390e82aaf1e1f0d3eade38cf81cce5809452..948129a68f47f73dd63f05a96dd11aee3a39ac2a 100644 --- a/src/common/micro/modd_param_ice.F90 +++ b/src/common/micro/modd_param_ice.F90 @@ -79,6 +79,8 @@ LOGICAL :: LSEDIM_AFTER ! sedimentation done before (.FALSE.) or after (.TRUE.) ! REAL :: XSPLIT_MAXCFL ! Maximum CFL number allowed for SPLIT scheme LOGICAL :: LSNOW_T ! Snow parameterization from Wurtz (2021) +! +LOGICAL :: LPACK_INTERP !To pack arrays before computing the different interpolations (kernels and other) END TYPE PARAM_ICE_t ! TYPE(PARAM_ICE_t), SAVE, TARGET :: PARAM_ICE @@ -98,7 +100,8 @@ LOGICAL, POINTER :: LWARM => NULL(), & LADJ_BEFORE => NULL(), & LADJ_AFTER => NULL(), & LSEDIM_AFTER => NULL(),& - LSNOW_T => NULL() + LSNOW_T => NULL(),& + LPACK_INTERP => NULL() REAL, POINTER :: XVDEPOSC => NULL(), & XFRACM90 => NULL(), & @@ -138,6 +141,7 @@ SUBROUTINE PARAM_ICE_ASSOCIATE() LADJ_AFTER => PARAM_ICE%LADJ_AFTER LSEDIM_AFTER => PARAM_ICE%LSEDIM_AFTER LSNOW_T => PARAM_ICE%LSNOW_T + LPACK_INTERP => PARAM_ICE%LPACK_INTERP ! XVDEPOSC => PARAM_ICE%XVDEPOSC XFRACM90 => PARAM_ICE%XFRACM90 diff --git a/src/common/micro/mode_ice4_fast_rg.F90 b/src/common/micro/mode_ice4_fast_rg.F90 index a3457b89e99265cb6324ac53302a3616e871d308..e861b521e69116ceee0607e9998e1debd2314a0e 100644 --- a/src/common/micro/mode_ice4_fast_rg.F90 +++ b/src/common/micro/mode_ice4_fast_rg.F90 @@ -97,15 +97,14 @@ REAL, DIMENSION(KPROMA, 8), INTENT(INOUT) :: PRG_TEND ! Individual tendencies INTEGER, PARAMETER :: IRCDRYG=1, IRIDRYG=2, IRIWETG=3, IRSDRYG=4, IRSWETG=5, IRRDRYG=6, & & IFREEZ1=7, IFREEZ2=8 LOGICAL, DIMENSION(KPROMA) :: GDRY, LLDRYG -INTEGER, DIMENSION(KPROMA) :: I1 INTEGER :: IGDRY -REAL, DIMENSION(KPROMA) :: ZVEC1, ZVEC2, ZVEC3 -INTEGER, DIMENSION(KPROMA) :: IVEC1, IVEC2 +REAL, DIMENSION(KPROMA) :: ZBUF1, ZBUF2, ZBUF3 +INTEGER, DIMENSION(KPROMA) :: IBUF1, IBUF2, IBUF3 REAL, DIMENSION(KPROMA) :: ZZW, & ZRDRYG_INIT, & !Initial dry growth rate of the graupeln ZRWETG_INIT !Initial wet growth rate of the graupeln REAL :: ZZW0D -INTEGER :: JJ, JL +INTEGER :: JL REAL(KIND=JPRB) :: ZHOOK_HANDLE !------------------------------------------------------------------------------- @@ -171,58 +170,23 @@ DO JL=1, KSIZE ENDDO ! Wet and dry collection of rs on graupel (6.2.1) -IGDRY = 0 DO JL = 1, KSIZE IF (PRST(JL)>ICED%XRTMIN(5) .AND. PRGT(JL)>ICED%XRTMIN(6) .AND. LDCOMPUTE(JL)) THEN - IGDRY = IGDRY + 1 - I1(IGDRY) = JL GDRY(JL) = .TRUE. ELSE GDRY(JL) = .FALSE. PRG_TEND(JL, IRSDRYG)=0. PRG_TEND(JL, IRSWETG)=0. - END IF + ENDIF ENDDO + IF(.NOT. LDSOFT) THEN + CALL INTERP_MICRO_2D(KPROMA, KSIZE, PLBDAG(:), PLBDAS(:), ICEP%NDRYLBDAG, ICEP%NDRYLBDAS, & + &ICEP%XDRYINTP1G, ICEP%XDRYINTP2G, ICEP%XDRYINTP1S, ICEP%XDRYINTP2S, & + &PARAMI%LPACK_INTERP, GDRY(:), IBUF1(:), IBUF2(:), IBUF3(:), ZBUF1(:), ZBUF2(:), ZBUF3(:), & + &IGDRY, & + &ICEP%XKER_SDRYG(:,:), ZZW(:)) IF(IGDRY>0)THEN - ! - !* 6.2.3 select the (PLBDAG,PLBDAS) couplet - ! - DO JJ = 1, IGDRY - ZVEC1(JJ) = PLBDAG(I1(JJ)) - ZVEC2(JJ) = PLBDAS(I1(JJ)) - END DO - ! - !* 6.2.4 find the next lower indice for the PLBDAG and for the PLBDAS - ! in the geometrical set of (Lbda_g,Lbda_s) couplet use to - ! tabulate the SDRYG-kernel - ! - ZVEC1(1:IGDRY)=MAX(1.00001, MIN(REAL(ICEP%NDRYLBDAG)-0.00001, & - ICEP%XDRYINTP1G*LOG(ZVEC1(1:IGDRY))+ICEP%XDRYINTP2G)) - IVEC1(1:IGDRY)=INT(ZVEC1(1:IGDRY) ) - ZVEC1(1:IGDRY)=ZVEC1(1:IGDRY)-REAL(IVEC1(1:IGDRY)) - ! - ZVEC2(1:IGDRY)=MAX(1.00001, MIN( REAL(ICEP%NDRYLBDAS)-0.00001, & - ICEP%XDRYINTP1S*LOG(ZVEC2(1:IGDRY))+ICEP%XDRYINTP2S)) - IVEC2(1:IGDRY)=INT(ZVEC2(1:IGDRY)) - ZVEC2(1:IGDRY)=ZVEC2(1:IGDRY)-REAL(IVEC2(1:IGDRY)) - ! - !* 6.2.5 perform the bilinear interpolation of the normalized - ! SDRYG-kernel - ! - DO JJ=1, IGDRY - ZVEC3(JJ) = ( ICEP%XKER_SDRYG(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_SDRYG(IVEC1(JJ)+1,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * ZVEC1(JJ) & - - ( ICEP%XKER_SDRYG(IVEC1(JJ) ,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_SDRYG(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - *(ZVEC1(JJ) - 1.0) - END DO - ZZW(:) = 0. - DO JJ = 1, IGDRY - ZZW(I1(JJ)) = ZVEC3(JJ) - END DO - ! !$mnh_expand_where(JL=1:KSIZE) WHERE(GDRY(1:KSIZE)) PRG_TEND(1:KSIZE, IRSWETG)=ICEP%XFSDRYG*ZZW(1:KSIZE) & ! RSDRYG @@ -245,11 +209,8 @@ ENDIF ! !* 6.2.6 accretion of raindrops on the graupeln ! -IGDRY = 0 DO JL = 1, KSIZE IF (PRRT(JL)>ICED%XRTMIN(3) .AND. PRGT(JL)>ICED%XRTMIN(6) .AND. LDCOMPUTE(JL)) THEN - IGDRY = IGDRY + 1 - I1(IGDRY) = JL GDRY(JL) = .TRUE. ELSE GDRY(JL) = .FALSE. @@ -258,45 +219,12 @@ DO JL = 1, KSIZE ENDDO IF(.NOT. LDSOFT) THEN ! + CALL INTERP_MICRO_2D(KPROMA, KSIZE, PLBDAG(:), PLBDAR(:), ICEP%NDRYLBDAG, ICEP%NDRYLBDAR, & + &ICEP%XDRYINTP1G, ICEP%XDRYINTP2G, ICEP%XDRYINTP1R, ICEP%XDRYINTP2R, & + &PARAMI%LPACK_INTERP, GDRY(:), IBUF1(:), IBUF2(:), IBUF3(:), ZBUF1(:), ZBUF2(:), ZBUF3(:), & + &IGDRY, & + &ICEP%XKER_RDRYG(:,:), ZZW(:)) IF(IGDRY>0) THEN - ! - !* 6.2.8 select the (PLBDAG,PLBDAR) couplet - ! - DO JJ = 1, IGDRY - ZVEC1(JJ) = PLBDAG(I1(JJ)) - ZVEC2(JJ) = PLBDAR(I1(JJ)) - ENDDO - ! - !* 6.2.9 find the next lower indice for the PLBDAG and for the PLBDAR - ! in the geometrical set of (Lbda_g,Lbda_r) couplet use to - ! tabulate the RDRYG-kernel - ! - ZVEC1(1:IGDRY)=MAX(1.00001, MIN( REAL(ICEP%NDRYLBDAG)-0.00001, & - ICEP%XDRYINTP1G*LOG(ZVEC1(1:IGDRY))+ICEP%XDRYINTP2G)) - IVEC1(1:IGDRY)=INT(ZVEC1(1:IGDRY)) - ZVEC1(1:IGDRY)=ZVEC1(1:IGDRY)-REAL(IVEC1(1:IGDRY)) - ! - ZVEC2(1:IGDRY)=MAX(1.00001, MIN( REAL(ICEP%NDRYLBDAR)-0.00001, & - ICEP%XDRYINTP1R*LOG(ZVEC2(1:IGDRY))+ICEP%XDRYINTP2R)) - IVEC2(1:IGDRY)=INT(ZVEC2(1:IGDRY)) - ZVEC2(1:IGDRY)=ZVEC2(1:IGDRY)-REAL(IVEC2(1:IGDRY)) - ! - !* 6.2.10 perform the bilinear interpolation of the normalized - ! RDRYG-kernel - ! - DO JJ=1, IGDRY - ZVEC3(JJ)= ( ICEP%XKER_RDRYG(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_RDRYG(IVEC1(JJ)+1,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * ZVEC1(JJ) & - - ( ICEP%XKER_RDRYG(IVEC1(JJ) ,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_RDRYG(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - *(ZVEC1(JJ) - 1.0) - END DO - ZZW(:) = 0. - DO JJ = 1, IGDRY - ZZW(I1(JJ)) = ZVEC3(JJ) - END DO - ! !$mnh_expand_where(JL=1:KSIZE) WHERE(GDRY(1:KSIZE)) PRG_TEND(1:KSIZE, IRRDRYG) = ICEP%XFRDRYG*ZZW(1:KSIZE) & ! RRDRYG @@ -439,6 +367,11 @@ DO JL=1, KSIZE ENDDO ! IF (LHOOK) CALL DR_HOOK('ICE4_FAST_RG', 1, ZHOOK_HANDLE) - +! +CONTAINS +! +INCLUDE "interp_micro.func.h" +! END SUBROUTINE ICE4_FAST_RG +! END MODULE MODE_ICE4_FAST_RG diff --git a/src/common/micro/mode_ice4_fast_rh.F90 b/src/common/micro/mode_ice4_fast_rh.F90 index 3d13263d73a3d9d3da053943622cc70cc540113b..c35275a392405ce6ac6d661d415b0f4191b5549b 100644 --- a/src/common/micro/mode_ice4_fast_rh.F90 +++ b/src/common/micro/mode_ice4_fast_rh.F90 @@ -93,9 +93,8 @@ INTEGER, PARAMETER :: IRCWETH=1, IRRWETH=2, IRIDRYH=3, IRIWETH=4, IRSDRYH=5, IRS & IFREEZ1=9, IFREEZ2=10 LOGICAL, DIMENSION(KPROMA) :: GWET INTEGER :: IGWET -INTEGER, DIMENSION(KPROMA) :: I1 -REAL, DIMENSION(KPROMA) :: ZVEC1, ZVEC2, ZVEC3 -INTEGER, DIMENSION(KPROMA) :: IVEC1, IVEC2 +REAL, DIMENSION(KPROMA) :: ZBUF1, ZBUF2, ZBUF3 +INTEGER, DIMENSION(KPROMA) :: IBUF1, IBUF2, IBUF3 REAL, DIMENSION(KPROMA) :: ZZW, & ZRDRYH_INIT, ZRWETH_INIT, & ZRDRYHG @@ -134,11 +133,8 @@ ENDDO ! !* 7.2.1 accretion of aggregates on the hailstones ! -IGWET = 0 DO JL = 1, KSIZE IF (PRHT(JL)>ICED%XRTMIN(7) .AND. PRST(JL)>ICED%XRTMIN(5) .AND. LDCOMPUTE(JL)) THEN - IGWET = IGWET + 1 - I1(IGWET) = JL GWET(JL) = .TRUE. ELSE GWET(JL) = .FALSE. @@ -147,45 +143,12 @@ DO JL = 1, KSIZE ENDIF ENDDO IF(.NOT. LDSOFT) THEN + CALL INTERP_MICRO_2D(KPROMA, KSIZE, PLBDAH(:), PLBDAS(:), ICEP%NWETLBDAH, ICEP%NWETLBDAS, & + &ICEP%XWETINTP1H, ICEP%XWETINTP2H, ICEP%XWETINTP1S, ICEP%XWETINTP2S, & + &PARAMI%LPACK_INTERP, GWET(:), IBUF1(:), IBUF2(:), IBUF3(:), ZBUF1(:), ZBUF2(:), ZBUF3(:), & + &IGWET, & + &ICEP%XKER_SWETH(:,:), ZZW(:)) IF(IGWET>0)THEN - ! - !* 7.2.3 select the (PLBDAH,PLBDAS) couplet - ! - DO JJ = 1, IGWET - ZVEC1(JJ) = PLBDAH(I1(JJ)) - ZVEC2(JJ) = PLBDAS(I1(JJ)) - ENDDO - ! - !* 7.2.4 find the next lower indice for the PLBDAG and for the PLBDAS - ! in the geometrical set of (Lbda_h,Lbda_s) couplet use to - ! tabulate the SWETH-kernel - ! - ZVEC1(1:IGWET) = MAX( 1.00001, MIN( REAL(ICEP%NWETLBDAH)-0.00001, & - ICEP%XWETINTP1H * LOG( ZVEC1(1:IGWET) ) + ICEP%XWETINTP2H ) ) - IVEC1(1:IGWET) = INT( ZVEC1(1:IGWET) ) - ZVEC1(1:IGWET) = ZVEC1(1:IGWET) - REAL( IVEC1(1:IGWET) ) - ! - ZVEC2(1:IGWET) = MAX( 1.00001, MIN( REAL(ICEP%NWETLBDAS)-0.00001, & - ICEP%XWETINTP1S * LOG( ZVEC2(1:IGWET) ) + ICEP%XWETINTP2S ) ) - IVEC2(1:IGWET) = INT( ZVEC2(1:IGWET) ) - ZVEC2(1:IGWET) = ZVEC2(1:IGWET) - REAL( IVEC2(1:IGWET) ) - ! - !* 7.2.5 perform the bilinear interpolation of the normalized - ! SWETH-kernel - ! - DO JJ = 1,IGWET - ZVEC3(JJ) = ( ICEP%XKER_SWETH(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_SWETH(IVEC1(JJ)+1,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * ZVEC1(JJ) & - - ( ICEP%XKER_SWETH(IVEC1(JJ) ,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_SWETH(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * (ZVEC1(JJ) - 1.0) - END DO - ZZW(:) = 0. - DO JJ = 1, IGWET - ZZW(I1(JJ)) = ZVEC3(JJ) - END DO - ! !$mnh_expand_where(JL=1:KSIZE) WHERE(GWET(1:KSIZE)) PRH_TEND(1:KSIZE, IRSWETH)=ICEP%XFSWETH*ZZW(1:KSIZE) & ! RSWETH @@ -207,11 +170,8 @@ ENDIF ! !* 7.2.6 accretion of graupeln on the hailstones ! -IGWET = 0 DO JL = 1, KSIZE IF (PRHT(JL)>ICED%XRTMIN(7) .AND. PRGT(JL)>ICED%XRTMIN(6) .AND. LDCOMPUTE(JL)) THEN - IGWET = IGWET + 1 - I1(IGWET) = JL GWET(JL) = .TRUE. ELSE GWET(JL) = .FALSE. @@ -220,45 +180,12 @@ DO JL = 1, KSIZE END IF ENDDO IF(.NOT. LDSOFT) THEN + CALL INTERP_MICRO_2D(KPROMA, KSIZE, PLBDAH(:), PLBDAG(:), ICEP%NWETLBDAH, ICEP%NWETLBDAG, & + &ICEP%XWETINTP1H, ICEP%XWETINTP2H, ICEP%XWETINTP1G, ICEP%XWETINTP2G, & + &PARAMI%LPACK_INTERP, GWET(:), IBUF1(:), IBUF2(:), IBUF3(:), ZBUF1(:), ZBUF2(:), ZBUF3(:), & + &IGWET, & + &ICEP%XKER_GWETH(:,:), ZZW(:)) IF(IGWET>0)THEN - ! - !* 7.2.8 select the (PLBDAH,PLBDAG) couplet - ! - DO JJ = 1, IGWET - ZVEC1(JJ) = PLBDAH(I1(JJ)) - ZVEC2(JJ) = PLBDAG(I1(JJ)) - END DO - ! - !* 7.2.9 find the next lower indice for the PLBDAH and for the PLBDAG - ! in the geometrical set of (Lbda_h,Lbda_g) couplet use to - ! tabulate the GWETH-kernel - ! - ZVEC1(1:IGWET) = MAX( 1.00001, MIN( REAL(ICEP%NWETLBDAG)-0.00001, & - ICEP%XWETINTP1H * LOG( ZVEC1(1:IGWET) ) + ICEP%XWETINTP2H ) ) - IVEC1(1:IGWET) = INT( ZVEC1(1:IGWET) ) - ZVEC1(1:IGWET) = ZVEC1(1:IGWET) - REAL( IVEC1(1:IGWET) ) - ! - ZVEC2(1:IGWET) = MAX( 1.00001, MIN( REAL(ICEP%NWETLBDAG)-0.00001, & - ICEP%XWETINTP1G * LOG( ZVEC2(1:IGWET) ) + ICEP%XWETINTP2G ) ) - IVEC2(1:IGWET) = INT( ZVEC2(1:IGWET) ) - ZVEC2(1:IGWET) = ZVEC2(1:IGWET) - REAL( IVEC2(1:IGWET) ) - ! - !* 7.2.10 perform the bilinear interpolation of the normalized - ! GWETH-kernel - ! - DO JJ = 1,IGWET - ZVEC3(JJ) = ( ICEP%XKER_GWETH(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_GWETH(IVEC1(JJ)+1,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * ZVEC1(JJ) & - - ( ICEP%XKER_GWETH(IVEC1(JJ) ,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_GWETH(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * (ZVEC1(JJ) - 1.0) - END DO - ZZW(:) = 0. - DO JJ = 1, IGWET - ZZW(I1(JJ)) = ZVEC3(JJ) - END DO - ! !$mnh_expand_where(JL=1:KSIZE) WHERE(GWET(1:KSIZE)) PRH_TEND(1:KSIZE, IRGWETH)=ICEP%XFGWETH*ZZW(1:KSIZE) & ! RGWETH @@ -279,11 +206,8 @@ ENDIF ! !* 7.2.11 accretion of raindrops on the hailstones ! -IGWET = 0 DO JL = 1, KSIZE IF (PRHT(JL)>ICED%XRTMIN(7) .AND. PRRT(JL)>ICED%XRTMIN(3) .AND. LDCOMPUTE(JL)) THEN - IGWET = IGWET + 1 - I1(IGWET) = JL GWET(JL) = .TRUE. ELSE GWET(JL) = .FALSE. @@ -291,45 +215,12 @@ DO JL = 1, KSIZE ENDIF ENDDO IF(.NOT. LDSOFT) THEN + CALL INTERP_MICRO_2D(KPROMA, KSIZE, PLBDAH(:), PLBDAR(:), ICEP%NWETLBDAH, ICEP%NWETLBDAR, & + &ICEP%XWETINTP1H, ICEP%XWETINTP2H, ICEP%XWETINTP1R, ICEP%XWETINTP2R, & + &PARAMI%LPACK_INTERP, GWET(:), IBUF1(:), IBUF2(:), IBUF3(:), ZBUF1(:), ZBUF2(:), ZBUF3(:), & + &IGWET, & + &ICEP%XKER_RWETH(:,:), ZZW(:)) IF(IGWET>0)THEN - ! - !* 7.2.12 select the (PLBDAH,PLBDAR) couplet - ! - DO JJ = 1, IGWET - ZVEC1(JJ) = PLBDAH(I1(JJ)) - ZVEC2(JJ) = PLBDAR(I1(JJ)) - ENDDO - ! - !* 7.2.13 find the next lower indice for the PLBDAH and for the PLBDAR - ! in the geometrical set of (Lbda_h,Lbda_r) couplet use to - ! tabulate the RWETH-kernel - ! - ZVEC1(1:IGWET)=MAX(1.00001, MIN( REAL(ICEP%NWETLBDAH)-0.00001, & - ICEP%XWETINTP1H*LOG(ZVEC1(1:IGWET))+ICEP%XWETINTP2H)) - IVEC1(1:IGWET)=INT(ZVEC1(1:IGWET)) - ZVEC1(1:IGWET)=ZVEC1(1:IGWET)-REAL(IVEC1(1:IGWET)) - ! - ZVEC2(1:IGWET)=MAX(1.00001, MIN( REAL(ICEP%NWETLBDAR)-0.00001, & - ICEP%XWETINTP1R*LOG(ZVEC2(1:IGWET))+ICEP%XWETINTP2R)) - IVEC2(1:IGWET)=INT(ZVEC2(1:IGWET)) - ZVEC2(1:IGWET)=ZVEC2(1:IGWET)-REAL(IVEC2(1:IGWET)) - ! - !* 7.2.14 perform the bilinear interpolation of the normalized - ! RWETH-kernel - ! - DO JJ=1, IGWET - ZVEC3(JJ)= ( ICEP%XKER_RWETH(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_RWETH(IVEC1(JJ)+1,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * ZVEC1(JJ) & - - ( ICEP%XKER_RWETH(IVEC1(JJ) ,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_RWETH(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - *(ZVEC1(JJ) - 1.0) - END DO - ZZW(:) = 0. - DO JJ = 1, IGWET - ZZW(I1(JJ)) = ZVEC3(JJ) - END DO - ! !$mnh_expand_where(JL=1:KSIZE) WHERE(GWET(1:KSIZE)) PRH_TEND(1:KSIZE, IRRWETH) = ICEP%XFRWETH*ZZW(1:KSIZE) & ! RRWETH @@ -476,5 +367,9 @@ ENDDO ! IF (LHOOK) CALL DR_HOOK('ICE4_FAST_RH', 1, ZHOOK_HANDLE) ! +CONTAINS +! +INCLUDE "interp_micro.func.h" +! END SUBROUTINE ICE4_FAST_RH END MODULE MODE_ICE4_FAST_RH diff --git a/src/common/micro/mode_ice4_fast_rs.F90 b/src/common/micro/mode_ice4_fast_rs.F90 index 6655b2061604da1a8fc9f135b27728317ad0ba93..d79105d0c9f6cd0b2de3959b14a4bc558aa16288 100644 --- a/src/common/micro/mode_ice4_fast_rs.F90 +++ b/src/common/micro/mode_ice4_fast_rs.F90 @@ -85,11 +85,11 @@ INTEGER, PARAMETER :: IRCRIMS=1, IRCRIMSS=2, IRSRIMCG=3, IRRACCS=4, IRRACCSS=5, & IFREEZ1=7, IFREEZ2=8 LOGICAL, DIMENSION(KPROMA) :: GRIM, GACC INTEGER :: IGRIM, IGACC -INTEGER, DIMENSION(KPROMA) :: I1 -REAL, DIMENSION(KPROMA) :: ZVEC1, ZVEC2, ZVEC3 -INTEGER, DIMENSION(KPROMA) :: IVEC1, IVEC2 -REAL, DIMENSION(KPROMA) :: ZZW, ZZW2, ZZW6, ZFREEZ_RATE +INTEGER, DIMENSION(KPROMA) :: IBUF1, IBUF2, IBUF3 +REAL, DIMENSION(KPROMA) :: ZBUF1, ZBUF2, ZBUF3 +REAL, DIMENSION(KPROMA) :: ZZW, ZZW1, ZZW2, ZZW3, ZFREEZ_RATE INTEGER :: JJ, JL +REAL :: ZZW0D REAL(KIND=JPRB) :: ZHOOK_HANDLE !------------------------------------------------------------------------------- ! @@ -137,11 +137,13 @@ ENDDO ! !* 5.1 cloud droplet riming of the aggregates ! -IGRIM = 0 DO JL=1, KSIZE IF (PRCT(JL)>ICED%XRTMIN(2) .AND. PRST(JL)>ICED%XRTMIN(5) .AND. LDCOMPUTE(JL)) THEN - IGRIM = IGRIM + 1 - I1(IGRIM) = JL +#if defined(REPRO48) || defined(REPRO55) + ZZW(JL) = PLBDAS(JL) +#else + ZZW(JL) = (PLBDAS(JL)**ICED%XALPHAS + ICED%XFVELOS**ICED%XALPHAS)**(1./ICED%XALPHAS) +#endif GRIM(JL) = .TRUE. ELSE GRIM(JL) = .FALSE. @@ -153,77 +155,29 @@ ENDDO ! ! Collection of cloud droplets by snow: this rate is used for riming (T<0) and for conversion/melting (T>0) IF(.NOT. LDSOFT) THEN + CALL INTERP_MICRO_1D(KPROMA, KSIZE, ZZW, ICEP%NGAMINC, ICEP%XRIMINTP1, ICEP%XRIMINTP2, & + PARAMI%LPACK_INTERP, GRIM(:), IBUF1, IBUF2, ZBUF1, ZBUF2, & + IGRIM, & + ICEP%XGAMINC_RIM1(:), ZZW1(:), ICEP%XGAMINC_RIM2(:), ZZW2(:), ICEP%XGAMINC_RIM4(:), ZZW3(:)) IF(IGRIM>0) THEN - ! - ! 5.1.1 select the PLBDAS - ! - DO JJ = 1, IGRIM -#if defined(REPRO48) || defined(REPRO55) - ZVEC1(JJ) = PLBDAS(I1(JJ)) -#else - ZVEC1(JJ) = (PLBDAS(I1(JJ))**ICED%XALPHAS + ICED%XFVELOS**ICED%XALPHAS)**(1./ICED%XALPHAS) -#endif - END DO - ! - ! 5.1.2 find the next lower indice for the PLBDAS in the geometrical - ! set of Lbda_s used to tabulate some moments of the incomplete - ! gamma function - ! - !$mnh_expand_where(JJ=1:IGRIM) - ZVEC2(1:IGRIM) = MAX( 1.00001, MIN( REAL(ICEP%NGAMINC)-0.00001, & - ICEP%XRIMINTP1 * LOG( ZVEC1(1:IGRIM) ) + ICEP%XRIMINTP2 ) ) - IVEC2(1:IGRIM) = INT( ZVEC2(1:IGRIM) ) - ZVEC2(1:IGRIM) = ZVEC2(1:IGRIM) - REAL( IVEC2(1:IGRIM) ) - ! - ! 5.1.3 perform the linear interpolation of the normalized - ! "2+XDS"-moment of the incomplete gamma function - ! - ZVEC1(1:IGRIM) = ICEP%XGAMINC_RIM1( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM) & - - ICEP%XGAMINC_RIM1( IVEC2(1:IGRIM) )*(ZVEC2(1:IGRIM) - 1.0) - !$mnh_end_expand_where(JJ=1:IGRIM) - ZZW(:) = 0. - DO JJ = 1, IGRIM - ZZW(I1(JJ)) = ZVEC1(JJ) - END DO ! ! 5.1.4 riming of the small sized aggregates ! !$mnh_expand_where(JL=1:KSIZE) WHERE (GRIM(1:KSIZE)) - PRS_TEND(1:KSIZE, IRCRIMSS) = ICEP%XCRIMSS * ZZW(1:KSIZE) * PRCT(1:KSIZE) & ! RCRIMSS + PRS_TEND(1:KSIZE, IRCRIMSS) = ICEP%XCRIMSS * ZZW1(1:KSIZE) * PRCT(1:KSIZE) & ! RCRIMSS #if defined(REPRO48) || defined(REPRO55) - * PLBDAS(1:KSIZE)**ICEP%XEXCRIMSS & + * PLBDAS(1:KSIZE)**ICEP%XEXCRIMSS & * PRHODREF(1:KSIZE)**(-ICED%XCEXVT) #else - * PRST(1:KSIZE)*(1+(ICED%XFVELOS/PLBDAS(1:KSIZE))**ICED%XALPHAS)**(-ICED%XNUS+ICEP%XEXCRIMSS/ICED%XALPHAS) & + * PRST(1:KSIZE)*(1+(ICED%XFVELOS/PLBDAS(1:KSIZE))**ICED%XALPHAS) & + **(-ICED%XNUS+ICEP%XEXCRIMSS/ICED%XALPHAS) & * PRHODREF(1:KSIZE)**(-ICED%XCEXVT+1.) & - * (PLBDAS(1:KSIZE)) ** (ICEP%XEXCRIMSS+ICED%XBS) + * (PLBDAS(1:KSIZE)) ** (ICEP%XEXCRIMSS+ICED%XBS) #endif END WHERE !$mnh_end_expand_where(JL=1:KSIZE) ! - ! 5.1.5 perform the linear interpolation of the normalized - ! "XBS"-moment of the incomplete gamma function (XGAMINC_RIM2) and - ! "XBG"-moment of the incomplete gamma function (XGAMINC_RIM4) - ! - !$mnh_expand_where(JJ=1:IGRIM) - ZVEC1(1:IGRIM) = ICEP%XGAMINC_RIM2( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM) & - - ICEP%XGAMINC_RIM2( IVEC2(1:IGRIM) )*(ZVEC2(1:IGRIM) - 1.0) - !$mnh_end_expand_where(JJ=1:IGRIM) - ZZW(:) = 0. - DO JJ = 1, IGRIM - ZZW(I1(JJ)) = ZVEC1(JJ) - END DO - - !$mnh_expand_where(JJ=1:IGRIM) - ZVEC1(1:IGRIM) = ICEP%XGAMINC_RIM4( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM) & - - ICEP%XGAMINC_RIM4( IVEC2(1:IGRIM) )*(ZVEC2(1:IGRIM) - 1.0) - !$mnh_end_expand_where(JJ=1:IGRIM) - ZZW2(:) = 0. - DO JJ = 1, IGRIM - ZZW2(I1(JJ)) = ZVEC1(JJ) - END DO - ! ! 5.1.6 riming-conversion of the large sized aggregates into graupeln ! ! @@ -234,7 +188,8 @@ IF(.NOT. LDSOFT) THEN * PLBDAS(1:KSIZE)**ICEP%XEXCRIMSG & * PRHODREF(1:KSIZE)**(-ICED%XCEXVT) #else - * PRST(1:KSIZE)*(1+(ICED%XFVELOS/PLBDAS(1:KSIZE))**(ICED%XALPHAS))**(-ICED%XNUS+ICEP%XEXCRIMSG/ICED%XALPHAS) & + * PRST(1:KSIZE)*(1+(ICED%XFVELOS/PLBDAS(1:KSIZE))**(ICED%XALPHAS)) & + **(-ICED%XNUS+ICEP%XEXCRIMSG/ICED%XALPHAS) & * PRHODREF(1:KSIZE)**(-ICED%XCEXVT+1.) & * PLBDAS(1:KSIZE)**(ICED%XBS+ICEP%XEXCRIMSG) #endif @@ -245,18 +200,20 @@ IF(.NOT. LDSOFT) THEN !Murakami 1990 !$mnh_expand_where(JL=1:KSIZE) WHERE(GRIM(1:KSIZE)) - ZZW6(1:KSIZE) = PRS_TEND(1:KSIZE, IRCRIMS) - PRS_TEND(1:KSIZE, IRCRIMSS) ! RCRIMSG + ZZW(1:KSIZE) = PRS_TEND(1:KSIZE, IRCRIMS) - PRS_TEND(1:KSIZE, IRCRIMSS) ! RCRIMSG #if defined(REPRO48) || defined(REPRO55) - PRS_TEND(1:KSIZE, IRSRIMCG)=ICEP%XSRIMCG * PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG*(1.0-ZZW(1:KSIZE)) + PRS_TEND(1:KSIZE, IRSRIMCG)=ICEP%XSRIMCG * PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG*(1.0-ZZW2(1:KSIZE)) #else - PRS_TEND(1:KSIZE, IRSRIMCG)=ICEP%XSRIMCG * PRST(1:KSIZE)*PRHODREF(1:KSIZE)*PLBDAS(1:KSIZE)**(ICEP%XEXSRIMCG+ICED%XBS)*(1.0-ZZW(1:KSIZE)) + PRS_TEND(1:KSIZE, IRSRIMCG)=ICEP%XSRIMCG * PRST(1:KSIZE)*PRHODREF(1:KSIZE) & + * PLBDAS(1:KSIZE)**(ICEP%XEXSRIMCG+ICED%XBS)*(1.0-ZZW2(1:KSIZE)) #endif - PRS_TEND(1:KSIZE, IRSRIMCG)=ZZW6(1:KSIZE)*PRS_TEND(1:KSIZE, IRSRIMCG)/ & + PRS_TEND(1:KSIZE, IRSRIMCG)=ZZW(1:KSIZE)*PRS_TEND(1:KSIZE, IRSRIMCG)/ & MAX(1.E-20, & #if defined(REPRO48) || defined(REPRO55) - ICEP%XSRIMCG3*ICEP%XSRIMCG2*PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG2*(1.-ZZW2(1:KSIZE)) - & + ICEP%XSRIMCG3*ICEP%XSRIMCG2*PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG2*(1.-ZZW3(1:KSIZE)) - & #else - ICEP%XSRIMCG3*ICEP%XSRIMCG2*PRST(1:KSIZE)*PRHODREF(1:KSIZE)*PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG2*(1.-ZZW2(1:KSIZE)) - & + ICEP%XSRIMCG3*ICEP%XSRIMCG2*PRST(1:KSIZE)*PRHODREF(1:KSIZE) & + *PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG2*(1.-ZZW3(1:KSIZE)) - & #endif ICEP%XSRIMCG3*PRS_TEND(1:KSIZE, IRSRIMCG)) END WHERE @@ -272,10 +229,10 @@ DO JL=1, KSIZE IF(GRIM(JL) .AND. PT(JL)<CST%XTT) THEN PRCRIMSS(JL)=MIN(ZFREEZ_RATE(JL), PRS_TEND(JL, IRCRIMSS)) ZFREEZ_RATE(JL)=MAX(0., ZFREEZ_RATE(JL)-PRCRIMSS(JL)) - ZZW(JL) = MIN(1., ZFREEZ_RATE(JL) / MAX(1.E-20, PRS_TEND(JL, IRCRIMS) - PRCRIMSS(JL))) ! proportion we are able to freeze - PRCRIMSG(JL) = ZZW(JL) * MAX(0., PRS_TEND(JL, IRCRIMS) - PRCRIMSS(JL)) ! RCRIMSG + ZZW0D = MIN(1., ZFREEZ_RATE(JL) / MAX(1.E-20, PRS_TEND(JL, IRCRIMS) - PRCRIMSS(JL))) ! proportion we are able to freeze + PRCRIMSG(JL) = ZZW0D * MAX(0., PRS_TEND(JL, IRCRIMS) - PRCRIMSS(JL)) ! RCRIMSG ZFREEZ_RATE(JL)=MAX(0., ZFREEZ_RATE(JL)-PRCRIMSG(JL)) - PRSRIMCG(JL) = ZZW(JL) * PRS_TEND(JL, IRSRIMCG) + PRSRIMCG(JL) = ZZW0D * PRS_TEND(JL, IRSRIMCG) PRSRIMCG(JL) = PRSRIMCG(JL) * MAX(0., -SIGN(1., -PRCRIMSG(JL))) PRCRIMSG(JL)=MAX(0., PRCRIMSG(JL)) @@ -288,11 +245,8 @@ ENDDO ! !* 5.2 rain accretion onto the aggregates ! -IGACC = 0 DO JL = 1, KSIZE IF (PRRT(JL)>ICED%XRTMIN(3) .AND. PRST(JL)>ICED%XRTMIN(5) .AND. LDCOMPUTE(JL)) THEN - IGACC = IGACC + 1 - I1(IGACC) = JL GACC(JL) = .TRUE. ELSE GACC(JL) = .FALSE. @@ -305,53 +259,17 @@ IF(.NOT. LDSOFT) THEN PRS_TEND(:, IRRACCS)=0. PRS_TEND(:, IRRACCSS)=0. PRS_TEND(:, IRSACCRG)=0. + CALL INTERP_MICRO_2D(KPROMA, KSIZE, PLBDAS, PLBDAR, ICEP%NACCLBDAS, ICEP%NACCLBDAR, & + &ICEP%XACCINTP1S, ICEP%XACCINTP2S, ICEP%XACCINTP1R, ICEP%XACCINTP2R,& + &PARAMI%LPACK_INTERP, GACC(:), IBUF1(:), IBUF2(:), IBUF3(:), ZBUF1(:), ZBUF2(:), ZBUF3(:), & + &IGACC, & + &ICEP%XKER_RACCSS(:,:), ZZW1(:), ICEP%XKER_RACCS(:,:), ZZW2(:), ICEP%XKER_SACCRG(:,:), ZZW3(:)) IF(IGACC>0)THEN - ! - ! - ! 5.2.1 select the (PLBDAS,PLBDAR) couplet - ! - DO JJ = 1, IGACC - ZVEC1(JJ) = PLBDAS(I1(JJ)) - ZVEC2(JJ) = PLBDAR(I1(JJ)) - ENDDO - ! - ! 5.2.2 find the next lower indice for the PLBDAS and for the PLBDAR - ! in the geometrical set of (Lbda_s,Lbda_r) couplet use to - ! tabulate the RACCSS-kernel - ! - !$mnh_expand_where(JJ=1:IGACC) - ZVEC1(1:IGACC) = MAX( 1.00001, MIN( REAL(ICEP%NACCLBDAS)-0.00001, & - ICEP%XACCINTP1S * LOG( ZVEC1(1:IGACC) ) + ICEP%XACCINTP2S ) ) - IVEC1(1:IGACC) = INT( ZVEC1(1:IGACC) ) - ZVEC1(1:IGACC) = ZVEC1(1:IGACC) - REAL( IVEC1(1:IGACC) ) - ! - ZVEC2(1:IGACC) = MAX( 1.00001, MIN( REAL(ICEP%NACCLBDAR)-0.00001, & - ICEP%XACCINTP1R * LOG( ZVEC2(1:IGACC) ) + ICEP%XACCINTP2R ) ) - IVEC2(1:IGACC) = INT( ZVEC2(1:IGACC) ) - ZVEC2(1:IGACC) = ZVEC2(1:IGACC) - REAL( IVEC2(1:IGACC) ) - !$mnh_end_expand_where(JJ=1:IGACC) - ! - ! 5.2.3 perform the bilinear interpolation of the normalized - ! RACCSS-kernel - ! - DO JJ = 1, IGACC - ZVEC3(JJ) = ( ICEP%XKER_RACCSS(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_RACCSS(IVEC1(JJ)+1,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * ZVEC1(JJ) & - - ( ICEP%XKER_RACCSS(IVEC1(JJ) ,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_RACCSS(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * (ZVEC1(JJ) - 1.0) - END DO - ZZW(:) = 0. - DO JJ = 1, IGACC - ZZW(I1(JJ)) = ZVEC3(JJ) - END DO - ! ! 5.2.4 raindrop accretion on the small sized aggregates ! !$mnh_expand_where(JL=1:KSIZE) WHERE(GACC(1:KSIZE)) - ZZW6(1:KSIZE) = & !! coef of RRACCS + ZZW(1:KSIZE) = & !! coef of RRACCS #if defined(REPRO48) || defined(REPRO55) ICEP%XFRACCSS*( PLBDAS(1:KSIZE)**ICED%XCXS )*( PRHODREF(1:KSIZE)**(-ICED%XCEXVT-1.) ) & #else @@ -360,52 +278,22 @@ IF(.NOT. LDSOFT) THEN *( ICEP%XLBRACCS1/((PLBDAS(1:KSIZE)**2) ) + & ICEP%XLBRACCS2/( PLBDAS(1:KSIZE) * PLBDAR(1:KSIZE) ) + & ICEP%XLBRACCS3/( (PLBDAR(1:KSIZE)**2)) )/PLBDAR(1:KSIZE)**4 - PRS_TEND(1:KSIZE, IRRACCSS) =ZZW(1:KSIZE)*ZZW6(1:KSIZE) + PRS_TEND(1:KSIZE, IRRACCSS) =ZZW1(1:KSIZE)*ZZW(1:KSIZE) END WHERE !$mnh_end_expand_where(JL=1:KSIZE) ! - ! 5.2.4b perform the bilinear interpolation of the normalized - ! RACCS-kernel - ! - DO JJ = 1, IGACC - ZVEC3(JJ) = ( ICEP%XKER_RACCS(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_RACCS(IVEC1(JJ)+1,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * ZVEC1(JJ) & - - ( ICEP%XKER_RACCS(IVEC1(JJ) ,IVEC2(JJ)+1)* ZVEC2(JJ) & - - ICEP%XKER_RACCS(IVEC1(JJ) ,IVEC2(JJ) )*(ZVEC2(JJ) - 1.0) ) & - * (ZVEC1(JJ) - 1.0) - END DO - ZZW(:) = 0. - DO JJ = 1, IGACC - ZZW(I1(JJ)) = ZVEC3(JJ) - END DO !$mnh_expand_where(JL=1:KSIZE) WHERE(GACC(1:KSIZE)) - PRS_TEND(1:KSIZE, IRRACCS) = ZZW(1:KSIZE)*ZZW6(1:KSIZE) + PRS_TEND(1:KSIZE, IRRACCS) = ZZW2(1:KSIZE)*ZZW(1:KSIZE) END WHERE !$mnh_end_expand_where(JL=1:KSIZE) - ! 5.2.5 perform the bilinear interpolation of the normalized - ! SACCRG-kernel - ! - DO JJ = 1, IGACC - ZVEC3(JJ) = ( ICEP%XKER_SACCRG(IVEC2(JJ)+1,IVEC1(JJ)+1)* ZVEC1(JJ) & - - ICEP%XKER_SACCRG(IVEC2(JJ)+1,IVEC1(JJ) )*(ZVEC1(JJ) - 1.0) ) & - * ZVEC2(JJ) & - - ( ICEP%XKER_SACCRG(IVEC2(JJ) ,IVEC1(JJ)+1)* ZVEC1(JJ) & - - ICEP%XKER_SACCRG(IVEC2(JJ) ,IVEC1(JJ) )*(ZVEC1(JJ) - 1.0) ) & - * (ZVEC2(JJ) - 1.0) - END DO - ZZW(:) = 0. - DO JJ = 1, IGACC - ZZW(I1(JJ)) = ZVEC3(JJ) - END DO ! ! 5.2.6 raindrop accretion-conversion of the large sized aggregates ! into graupeln ! !$mnh_expand_where(JL=1:KSIZE) WHERE(GACC(1:KSIZE)) - PRS_TEND(1:KSIZE, IRSACCRG) = ICEP%XFSACCRG*ZZW(1:KSIZE)* & ! RSACCRG + PRS_TEND(1:KSIZE, IRSACCRG) = ICEP%XFSACCRG*ZZW3(1:KSIZE)* & ! RSACCRG #if defined(REPRO48) || defined(REPRO55) ( PLBDAS(1:KSIZE)**(ICED%XCXS-ICED%XBS) )*( PRHODREF(1:KSIZE)**(-ICED%XCEXVT-1.) ) & #else @@ -484,5 +372,9 @@ ENDDO IF (LHOOK) CALL DR_HOOK('ICE4_FAST_RS', 1, ZHOOK_HANDLE) ! +CONTAINS +! +INCLUDE "interp_micro.func.h" +! END SUBROUTINE ICE4_FAST_RS END MODULE MODE_ICE4_FAST_RS diff --git a/src/common/micro/mode_ice4_rsrimcg_old.F90 b/src/common/micro/mode_ice4_rsrimcg_old.F90 index 865ed74004ea51dcc1a28356f2046e03a0413705..970f214f322bc79eec5ad895d366595959e0974e 100644 --- a/src/common/micro/mode_ice4_rsrimcg_old.F90 +++ b/src/common/micro/mode_ice4_rsrimcg_old.F90 @@ -6,7 +6,7 @@ MODULE MODE_ICE4_RSRIMCG_OLD IMPLICIT NONE CONTAINS -SUBROUTINE ICE4_RSRIMCG_OLD(CST, ICEP, ICED, KPROMA, KSIZE, LDSOFT, LDCOMPUTE, & +SUBROUTINE ICE4_RSRIMCG_OLD(CST, PARAMI, ICEP, ICED, KPROMA, KSIZE, LDSOFT, LDCOMPUTE, & &PRHODREF, & &PLBDAS, & &PT, PRCT, PRST, & @@ -31,6 +31,7 @@ SUBROUTINE ICE4_RSRIMCG_OLD(CST, ICEP, ICED, KPROMA, KSIZE, LDSOFT, LDCOMPUTE, & ! ------------ ! USE MODD_CST, ONLY: CST_t +USE MODD_PARAM_ICE, ONLY: PARAM_ICE_t USE MODD_RAIN_ICE_DESCR, ONLY: RAIN_ICE_DESCR_t USE MODD_RAIN_ICE_PARAM, ONLY: RAIN_ICE_PARAM_t USE PARKIND1, ONLY : JPRB @@ -41,6 +42,7 @@ IMPLICIT NONE !* 0.1 Declarations of dummy arguments : ! TYPE(CST_t), INTENT(IN) :: CST +TYPE(PARAM_ICE_t), INTENT(IN) :: PARAMI TYPE(RAIN_ICE_PARAM_t), INTENT(IN) :: ICEP TYPE(RAIN_ICE_DESCR_t), INTENT(IN) :: ICED INTEGER, INTENT(IN) :: KPROMA, KSIZE @@ -57,8 +59,8 @@ REAL, DIMENSION(KPROMA), INTENT(OUT) :: PRSRIMCG_MR ! Mr change due to c ! LOGICAL, DIMENSION(KPROMA) :: GRIM INTEGER :: IGRIM -REAL, DIMENSION(KPROMA) :: ZVEC1, ZVEC2 -INTEGER, DIMENSION(KPROMA) :: IVEC1, IVEC2 +REAL, DIMENSION(KPROMA) :: ZBUF1, ZBUF2 +INTEGER, DIMENSION(KPROMA) :: IBUF1, IBUF2 REAL, DIMENSION(KPROMA) :: ZZW INTEGER :: JL REAL(KIND=JPRB) :: ZHOOK_HANDLE @@ -73,49 +75,15 @@ IF (LHOOK) CALL DR_HOOK('ICE4_RSRIMCG_OLD', 0, ZHOOK_HANDLE) PRSRIMCG_MR(:)=0. ! IF(.NOT. LDSOFT) THEN - IGRIM = 0 DO JL = 1, KSIZE - IF(PRCT(JL)>ICED%XRTMIN(2) .AND. PRST(JL)>ICED%XRTMIN(5) .AND. LDCOMPUTE(JL) .AND. PT(JL)<CST%XTT) THEN - IGRIM = IGRIM + 1 - IVEC1(IGRIM) = JL - GRIM(JL) = .TRUE. - ELSE - GRIM(JL) = .FALSE. - ENDIF + GRIM(JL)=PRCT(JL)>ICED%XRTMIN(2) .AND. PRST(JL)>ICED%XRTMIN(5) .AND. LDCOMPUTE(JL) .AND. PT(JL)<CST%XTT ENDDO + CALL INTERP_MICRO_1D(KPROMA, KSIZE, PLBDAS(:), ICEP%NGAMINC, ICEP%XRIMINTP1, ICEP%XRIMINTP2, & + &PARAMI%LPACK_INTERP, GRIM(:), IBUF1, IBUF2, ZBUF1, ZBUF2, & + &IGRIM, & + &ICEP%XGAMINC_RIM2, ZZW) ! IF(IGRIM>0) THEN - ! - ! 5.1.1 select the PLBDAS - ! - DO JL = 1, IGRIM - ZVEC1(JL) = PLBDAS(IVEC1(JL)) - ENDDO - ! - ! 5.1.2 find the next lower indice for the PLBDAS in the geometrical - ! set of Lbda_s used to tabulate some moments of the incomplete - ! gamma function - ! - ZVEC2(1:IGRIM) = MAX( 1.00001, MIN(REAL(ICEP%NGAMINC)-0.00001, & - ICEP%XRIMINTP1 * LOG( ZVEC1(1:IGRIM) ) + ICEP%XRIMINTP2 ) ) - IVEC2(1:IGRIM) = INT( ZVEC2(1:IGRIM) ) - ZVEC2(1:IGRIM) = ZVEC2(1:IGRIM) - REAL( IVEC2(1:IGRIM) ) - - ! - ! 5.1.5 perform the linear interpolation of the normalized - ! "XBS"-moment of the incomplete gamma function (XGAMINC_RIM2) - ! - ZVEC1(1:IGRIM) = ICEP%XGAMINC_RIM2( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM) & - - ICEP%XGAMINC_RIM2( IVEC2(1:IGRIM) )*(ZVEC2(1:IGRIM) - 1.0) - ZZW(:) = 0. - DO JL = 1, IGRIM - ZZW(IVEC1(JL)) = ZVEC1(JL) - ENDDO - - ! - ! 5.1.6 riming-conversion of the large sized aggregates into graupeln - ! - ! !$mnh_expand_where(JL=1:KSIZE) WHERE(GRIM(1:KSIZE)) PRSRIMCG_MR(1:KSIZE) = ICEP%XSRIMCG * PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG & ! RSRIMCG @@ -132,5 +100,9 @@ ENDIF ! IF (LHOOK) CALL DR_HOOK('ICE4_RSRIMCG_OLD', 1, ZHOOK_HANDLE) ! +CONTAINS +! +INCLUDE "interp_micro.func.h" +! END SUBROUTINE ICE4_RSRIMCG_OLD END MODULE MODE_ICE4_RSRIMCG_OLD diff --git a/src/common/micro/mode_ice4_tendencies.F90 b/src/common/micro/mode_ice4_tendencies.F90 index d6eebfdf81ace36fcae6b31c221eedf33cadf132..fd722114c661d755f28eb062e2075deac77b44c0 100644 --- a/src/common/micro/mode_ice4_tendencies.F90 +++ b/src/common/micro/mode_ice4_tendencies.F90 @@ -255,7 +255,7 @@ ELSE ZLBDAS(1:KSIZE)=0. END WHERE !$mnh_end_expand_where(JL=1:KSIZE) - CALL ICE4_RSRIMCG_OLD(CST, ICEP, ICED, KPROMA, KSIZE, ODSOFT, LDCOMPUTE, & + CALL ICE4_RSRIMCG_OLD(CST, PARAMI, ICEP, ICED, KPROMA, KSIZE, ODSOFT, LDCOMPUTE, & &PRHODREF, & &ZLBDAS, & &ZT, ZVART(:,IRC), ZVART(:,IRS), & diff --git a/src/mesonh/ext/default_desfmn.f90 b/src/mesonh/ext/default_desfmn.f90 index ade773022d23c3b0a5a10956b158429b0fde831d..66209eacac13278222d72a1c75636e1f653d9bc0 100644 --- a/src/mesonh/ext/default_desfmn.f90 +++ b/src/mesonh/ext/default_desfmn.f90 @@ -880,6 +880,7 @@ IF (KMI == 1) THEN LDEPOSC = .FALSE. XVDEPOSC= 0.02 ! 2 cm/s LSNOW_T=.FALSE. + LPACK_INTERP=.TRUE. END IF ! !------------------------------------------------------------------------------- diff --git a/src/testprogs/rain_ice/main_rain_ice.F90 b/src/testprogs/rain_ice/main_rain_ice.F90 index 2f1fa18699aec824bee618e5a4e41288f6e15572..12c5e951e1644d3abf807a99141ef0fb44d1272a 100644 --- a/src/testprogs/rain_ice/main_rain_ice.F90 +++ b/src/testprogs/rain_ice/main_rain_ice.F90 @@ -407,6 +407,7 @@ LSEDIM_AFTER=.FALSE. ! Sedimentation done after microphysics XSPLIT_MAXCFL=0.8 LDEPOSC=.FALSE. ! water deposition on vegetation XVDEPOSC=0.02 ! deposition speed (2 cm.s-1) +LPACK_INTERP=.TRUE. ! ! 2. Set implicit default values for MODD_RAIN_ICE_DESCR ! et MODD_RAIN_ICE_PARAM