Index: LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.F90 (revision 3197) +++ LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.F90 (revision 3198) @@ -20,4 +20,7 @@ use dimphy, only: klon, klev use assert_m, only: assert + USE ioipsl_getin_p_mod, ONLY : getin_p + USE vertical_layers_mod, ONLY : presnivs + include "YOMCST.h" include "clesphys.h" @@ -111,4 +114,27 @@ REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels REAL BVSEC ! Security to avoid negative BVF + + REAL, DIMENSION(klev+1) ::HREF + LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true. + LOGICAL, SAVE :: firstcall = .TRUE. + !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp) + + CHARACTER (LEN=20) :: modname='flott_gwd_rando' + CHARACTER (LEN=80) :: abort_message + + + + IF (firstcall) THEN + ! Cle introduite pour resoudre un probleme de non reproductibilite + ! Le but est de pouvoir tester de revenir a la version precedenete + ! A eliminer rapidement + CALL getin_p('gwd_reproductibilite_mpiomp',gwd_reproductibilite_mpiomp) + IF (NW+4*(NA-1)+NA>=KLEV) THEN + abort_message = 'NW+3*NA>=KLEV Probleme pour generation des ondes' + CALL abort_physic (modname,abort_message,1) + ENDIF + firstcall=.false. +! CALL iophys_ini + ENDIF !----------------------------------------------------------------- @@ -205,12 +231,26 @@ ! Launching altitude + IF (gwd_reproductibilite_mpiomp) THEN + ! Reprend la formule qui calcule PH en fonction de PP=play + DO LL = 2, KLEV + HREF(LL) = EXP((LOG(presnivs(LL)) + LOG(presnivs(LL - 1))) / 2.) + end DO + HREF(KLEV + 1) = 0. + HREF(1) = 2. * presnivs(1) - HREF(2) + ELSE + HREF(1:KLEV)=PH(KLON/2,1:KLEV) + ENDIF + LAUNCH=0 LTROP =0 DO LL = 1, KLEV - IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XLAUNCH) LAUNCH = LL + IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL ENDDO DO LL = 1, KLEV - IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XTROP) LTROP = LL + IF (HREF(LL) / HREF(1) > XTROP) LTROP = LL ENDDO + !LAUNCH=22 ; LTROP=33 +! print*,'LAUNCH=',LAUNCH,'LTROP=',LTROP + ! PRINT *,'LAUNCH IN ACAMARA:',LAUNCH @@ -293,8 +333,5 @@ JW = 0 - DO JP = 1, NP - DO JK = 1, NK - DO JO = 1, NO - JW = JW + 1 + DO JW = 1, NW ! Angle DO II = 1, KLON @@ -340,6 +377,4 @@ ! RUW0(JW, II) = RUWFRT ENDDO - end DO - end DO end DO Index: LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.F90 (revision 3197) +++ LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.F90 (revision 3198) @@ -18,4 +18,7 @@ use dimphy, only: klon, klev use assert_m, only: assert + USE ioipsl_getin_p_mod, ONLY : getin_p + USE vertical_layers_mod, ONLY : presnivs + include "YOMCST.h" include "clesphys.h" @@ -103,7 +106,31 @@ REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels REAL PSEC ! Security to avoid division by 0 pressure - REAL PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels REAL BVSEC ! Security to avoid negative BVF + REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3 + + REAL, DIMENSION(klev+1) ::HREF + + LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true. + LOGICAL, SAVE :: firstcall = .TRUE. + !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp) + + CHARACTER (LEN=20) :: modname='flott_gwd_rando' + CHARACTER (LEN=80) :: abort_message + + + + IF (firstcall) THEN + ! Cle introduite pour resoudre un probleme de non reproductibilite + ! Le but est de pouvoir tester de revenir a la version precedenete + ! A eliminer rapidement + CALL getin_p('gwd_reproductibilite_mpiomp',gwd_reproductibilite_mpiomp) + IF (NW+3*NA>=KLEV) THEN + abort_message = 'NW+3*NA>=KLEV Probleme pour generation des ondes' + CALL abort_physic (modname,abort_message,1) + ENDIF + firstcall=.false. + ENDIF + !----------------------------------------------------------------- @@ -156,4 +183,5 @@ ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ +IF (1==0) THEN !ONLINE call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), & @@ -167,4 +195,5 @@ "FLOTT_GWD_RANDO klev") !END ONLINE +ENDIF IF(DELTAT < DTIME)THEN @@ -183,21 +212,32 @@ DO LL = 2, KLEV PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.) - PHM1(:, LL) = 1. / PH(:, LL) - end DO - + end DO PH(:, KLEV + 1) = 0. - PHM1(:, KLEV + 1) = 1. / PSEC PH(:, 1) = 2. * PP(:, 1) - PH(:, 2) ! Launching altitude + + !Pour revenir a la version non reproductible en changeant le nombre de process + IF (gwd_reproductibilite_mpiomp) THEN + ! Reprend la formule qui calcule PH en fonction de PP=play + DO LL = 2, KLEV + HREF(LL) = EXP((LOG(presnivs(LL)) + LOG(presnivs(LL - 1))) / 2.) + end DO + HREF(KLEV + 1) = 0. + HREF(1) = 2. * presnivs(1) - HREF(2) + ELSE + HREF(1:KLEV)=PH(KLON/2,1:KLEV) + ENDIF LAUNCH=0 LTROP =0 DO LL = 1, KLEV - IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XLAUNCH) LAUNCH = LL + IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL ENDDO DO LL = 1, KLEV - IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XTROP) LTROP = LL + IF (HREF(LL) / HREF(1) > XTROP) LTROP = LL ENDDO + !LAUNCH=22 ; LTROP=33 +! print*,'LAUNCH=',LAUNCH,'LTROP=',LTROP ! Log pressure vert. coordinate @@ -245,21 +285,20 @@ ! waves characteristics in an almost stochastic way - JW = 0 - DO JP = 1, NP - DO JK = 1, NK - DO JO = 1, NO - JW = JW + 1 + DO JW = 1, NW ! Angle DO II = 1, KLON ! Angle (0 or PI so far) - ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) & + RAN_NUM_1=MOD(TT(II, JW) * 10., 1.) + RAN_NUM_2= MOD(TT(II, JW) * 100., 1.) + ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) & * RPI / 2. ! Horizontal wavenumber amplitude - ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.) + ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2 ! Horizontal phase speed CPHA = 0. DO JJ = 1, NA + RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.) CPHA = CPHA + & - CMAX*2.*(MOD(TT(II, JW+3*JJ)**2, 1.)-0.5)*SQRT(3.)/SQRT(NA*1.) + CMAX*2.*(RAN_NUM_3 -0.5)*SQRT(3.)/SQRT(NA*1.) END DO IF (CPHA.LT.0.) THEN @@ -276,7 +315,5 @@ RUW0(JW, II) = RUWMAX ENDDO - end DO - end DO - end DO + ENDDO ! 4. COMPUTE THE FLUXES @@ -417,4 +454,5 @@ ENDDO + END SUBROUTINE FLOTT_GWD_RANDO Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 3197) +++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 3198) @@ -1294,5 +1294,5 @@ ENDDO !!! jyg le 07/02/2012 et le 10/04/2013 - DO k = 1, klev + DO k = 1, klev+1 DO j = 1, knon i = ni(j) @@ -1300,5 +1300,10 @@ !! ytke(j,k) = tke(i,k,nsrf) ytke(j,k) = tke_x(i,k,nsrf) + ENDDO + ENDDO !>jyg + DO k = 1, klev + DO j = 1, knon + i = ni(j) !FC y_treedrg(j,k) = treedrg(i,k,nsrf) @@ -2408,5 +2413,5 @@ IF (iflag_split .eq.0) THEN wake_dltke(:,:,nsrf) = 0. - DO k = 1, klev + DO k = 1, klev+1 DO j = 1, knon i = ni(j) @@ -2421,5 +2426,5 @@ ELSE ! (iflag_split .eq.0) - DO k = 1, klev + DO k = 1, klev+1 DO j = 1, knon i = ni(j) Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 3197) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 3198) @@ -4390,5 +4390,5 @@ - CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pbl_tke) + CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke) Index: LMDZ6/trunk/libf/phylmd/tend_to_tke.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/tend_to_tke.F90 (revision 3197) +++ LMDZ6/trunk/libf/phylmd/tend_to_tke.F90 (revision 3198) @@ -32,5 +32,5 @@ !************************************************************************************** - SUBROUTINE tend_to_tke(dt,plev,exner,temp,windu,windv,dt_a,du_a,dv_a,tke) + SUBROUTINE tend_to_tke(dt,plev,exner,temp,windu,windv,dt_a,du_a,dv_a,pctsrf,tke) USE dimphy, ONLY: klon, klev @@ -53,4 +53,5 @@ REAL du_a(klon,klev) ! Zonal wind speed tendency [m/s], grid-cell average or for a one subsurface REAL dv_a(klon,klev) ! Meridional wind speed tendency [m/s], grid-cell average or for a one subsurface + REAL pctsrf(klon,nbsrf+1) ! Turbulent Kinetic energy [m2/s2], grid-cell average or for a subsurface ! Inputs/Outputs @@ -119,14 +120,14 @@ - DO isrf=1,nbsrf + DO isrf=1,nsrf DO k=1,klev - tke(:,k,isrf)= tke(:,k,isrf)+tendu(:,k)+tendv(:,k)+tendt(:,k) - tke(:,k,isrf)= max(tke(:,k,isrf),1.e-10) + DO i=1,klon + IF (pctsrf(i,isrf)>0.) THEN + tke(i,k,isrf)= tke(i,k,isrf)+tendu(i,k)+tendv(i,k)+tendt(i,k) + tke(i,k,isrf)= max(tke(i,k,isrf),1.e-10) + ENDIF + ENDDO ENDDO ENDDO - -! dtke_t(:,:)=tendt(:,:) -! dtke_u(:,:)=tendu(:,:) -! dtke_v(:,:)=tendv(:,:)