MODULE lmdz_wake_pkupper PUBLIC wake_pkupper CONTAINS SUBROUTINE wake_pkupper (klon, klev, ptop, ph, p, pupper, kupper, & dth, hw_, rho, delta_t_min_in, & ktop, wk_adv, h_zzz, ptop1, ktop1) USE lmdz_wake_ini , ONLY : wk_pupper USE lmdz_wake_ini , ONLY : RG USE lmdz_wake_ini , ONLY : hwmin USE lmdz_wake_ini , ONLY : iflag_wk_new_ptop, wk_delta_t_min, wk_frac_int_delta_t USE lmdz_wake_ini , ONLY : wk_int_delta_t_min IMPLICIT NONE INTEGER, INTENT(IN) :: klon,klev REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: ph, p REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: rho LOGICAL, DIMENSION (klon) , INTENT(IN) :: wk_adv REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: dth REAL, INTENT(IN) :: delta_t_min_in REAL, DIMENSION (klon) , INTENT(OUT) :: hw_ REAL, DIMENSION (klon) , INTENT(OUT) :: ptop INTEGER, DIMENSION (klon) , INTENT(OUT) :: Ktop REAL, DIMENSION (klon) , INTENT(OUT) :: pupper INTEGER, DIMENSION (klon) , INTENT(OUT) :: kupper REAL, DIMENSION (klon) , INTENT(OUT) :: h_zzz !! REAL, DIMENSION (klon) , INTENT(OUT) :: Ptop1 !! INTEGER, DIMENSION (klon) , INTENT(OUT) :: ktop1 !! INTEGER :: i,k LOGICAL, DIMENSION (klon) :: wk_active REAL :: delta_t_min REAL, DIMENSION (klon) :: dthmin REAL, DIMENSION (klon) :: ptop_provis,ptop_new REAL, DIMENSION (klon) :: z, dz REAL, DIMENSION (klon) :: sum_dth INTEGER, DIMENSION (klon) :: k_ptop_provis REAL, DIMENSION (klon) :: zk_ptop_provis REAL, DIMENSION (klon) :: omega !! REAL, DIMENSION (klon,klev+1) :: int_dth !! REAL, DIMENSION (klon,klev+1) :: dzz !! REAL, DIMENSION (klon,klev+1) :: zzz !! REAL, DIMENSION (klon) :: frac_int_dth !! REAL :: ddd!! REAL :: www INTEGER, SAVE :: ipas=0 !$OMP THREADPRIVATE(ipas) !INTEGER, SAVE :: compte=0 ! LJYF : a priori z, dz sum_dth sont aussi des variables internes ! Les eliminer apres verification convergence numerique !compte=compte+1 !print*,'compte=',compte ! Determine Ptop from buoyancy integral ! --------------------------------------- ! - 1/ Pressure of the level where dth changes sign. !print*,'WAKE LJYF' if (iflag_wk_new_ptop==0) then delta_t_min=delta_t_min_in else delta_t_min=wk_delta_t_min endif DO i = 1, klon ptop_provis(i) = ph(i, 1) k_ptop_provis(i) = 1 END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) k_ptop_provis(i) = k END IF END DO END DO ! - 2/ dth integral DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END IF END DO END DO ! - 3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon IF (wk_adv(i)) THEN hw_(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw_(i) = amax1(hwmin, hw_(i)) END IF END DO ! - 4/ now, get Ptop DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd ktop(i) = 0 z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw_(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i) ktop(i) = k END IF END IF END DO END DO ! 4.5/Correct ktop and ptop DO i = 1, klon ptop_new(i) = ptop(i) END DO DO k = klev, 2, -1 DO i = 1, klon ! IM v3JYG; IF (k .GE. ktop(i) IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon ptop(i) = ptop_new(i) END DO DO k = klev, 1, -1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (ph(i,k+1)<ptop(i)) ktop(i) = k END IF END DO END DO ! IF (prt_level>=10) THEN ! PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) ! ENDIF ! ----------------------------------------------------------------------- ! nouveau calcul de hw et ptop ! ----------------------------------------------------------------------- !if (iflag_wk_new_ptop>0) then do i=1,klon ptop1(i)=ph(i,1) ktop1(i)=1 h_zzz(i)=0. enddo IF (iflag_wk_new_ptop/=0) THEN int_dth(1:klon,1:klev+1)=0. DO i = 1, klon IF (wk_adv(i)) THEN int_dth(i,1) = 0. END IF END DO if (abs(iflag_wk_new_ptop) == 1 ) then DO k = 2, klev+1 Do i = 1, klon IF (wk_adv(i)) THEN if (k<=k_ptop_provis(i)) then ddd=dth(i,k-1)*(ph(i,k-1) - max(ptop_provis(i),ph(i,k))) !ddd=dth(i,k-1)*(ph(i,k-1) - ph(i,k)) else ddd=0. endif int_dth(i,k) = int_dth(i,k-1) + ddd !ELSE ! int_dth(i,k) = 0. END IF END DO END DO else k_ptop_provis(:)=klev+1 dthmin(:)=dth(:,1) ! calcul de l'int??grale de dT * dP jusqu'au dernier ! niveau avec dT<0. (en s'assurant qu'on a bien un ! dT negatif plus bas) DO k = 1, klev DO i = 1, klon dthmin(i)=min(dthmin(i),dth(i,k)) ddd=dth(i,k)*(ph(i,k)-ph(i,k+1)) if (dthmin(i)<0.) then if (k>=k_ptop_provis(i)) then ddd=0. else if (dth(i,k)>=0.) then ddd=0. k_ptop_provis(i)=k+1 endif endif int_dth(i,k+1) = int_dth(i,k)+ ddd ENDDO ENDDO DO i = 1, klon if ( k_ptop_provis(i)==klev+1 .or. .not. wk_adv(i)) then k_ptop_provis(i)=1 endif ENDDO endif ! (abs(iflag_wk_new_ptop) == 1 ) ! print*, 'xxx, int_dth', (k,int_dth(1,k),k=1,klev) ! print*, 'xxx, k_ptop_provis', k_ptop_provis(1) ! On se limite ?? des poches avec integrale dT * dp < -wk_int_delta_t_min do i=1,klon if (int_dth(i,k_ptop_provis(i)) > -wk_int_delta_t_min .or. k_ptop_provis(i)==1) then !if (1==0) then wk_active(i)=.false. ptop(i)=ph(i,1) ktop(i)=1 hw_(i)=0. else wk_active(i)=wk_adv(i) endif enddo DO i=1,klon IF (wk_active(i)) THEN frac_int_dth(i)=wk_frac_int_delta_t*int_dth(i,k_ptop_provis(i)) ENDIF ENDDO DO k = 1,klev DO i =1, klon IF (wk_active(i)) THEN IF (int_dth(i,k)>=frac_int_dth(i)) THEN ktop1(i) = min(k, k_ptop_provis(i)) !ktop1(i) = k !print*,ipas,'yyy ktop1= ',ktop1 ENDIF ENDIF END DO END DO !print*, 'LAMINE' DO i = 1, klon IF (wk_active(i)) THEN !print*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop1 ddd=int_dth(i,ktop1(i)+1)-int_dth(i,ktop1(i)) if (ddd==0.) then omega(i)=0. else omega(i) = (frac_int_dth(i) - int_dth(i,ktop1(i)))/ddd endif !! print*,'OMEGA ',omega(i) END IF END DO !! print*, 'xxx' DO i = 1, klon IF (wk_active(i)) THEN ! print*, 'xxx, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ', & ! int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ! print*, 'xxx, omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ', & !e omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ptop1(i) = min((1 - omega(i))*ph(i,ktop1(i)) + omega(i)*ph(i,ktop1(i)+1), ph(i,1)) END IF END DO DO i=1, klon IF (wk_active(i)) THEN zzz(i, 1) = 0 END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_active(i)) THEN dzz(i,k) = (ph(i,k) - ph(i,k+1))/(rho(i,k)*RG) zzz(i,k+1) = zzz(i,k) + dzz(i,k) END IF END DO END DO DO i =1, klon IF (wk_active(i)) THEN h_zzz(i) = max((1- omega(i))*zzz(i,ktop1(i)) + omega(i)*zzz(i,ktop1(i)+1), hwmin) END IF END DO ENDIF ! (iflag_wk_new_ptop/=0) !if (iflag_wk_new_ptop==2) then IF (iflag_wk_new_ptop>0) THEN do i=1,klon ptop(i)=ptop1(i) ktop(i)=ktop1(i) hw_(i)=h_zzz(i) enddo !endif ENDIF kupper = 0 IF (0.<wk_pupper .and. wk_pupper<1.) THEN ! Choose an integration bound well above wake top ! ----------------------------------------------------------------- ! Pupper = 50000. ! melting level ! Pupper = 60000. ! Pupper = 80000. ! essais pour case_e DO i = 1, klon ! pupper(i) = 0.6*ph(i, 1) pupper(i) = wk_pupper*ph(i, 1) pupper(i) = max(pupper(i), 45000.) ! cc Pupper(i) = 60000. END DO ELSE IF (1.<=wk_pupper) THEN DO i=1, klon ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) ! pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.) pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-5000.) END DO ELSE www=abs(wk_pupper) DO i=1, klon ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) ! pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.) pupper(i) = min( exp(www*log(ptop(i))+(1.-www)*log(ph(i, 1)) ) , ptop(i)-5000.) END DO END IF ! -5/ Determination de kupper DO k = klev, 1, -1 DO i = 1, klon IF (ph(i,k+1)<pupper(i)) kupper(i) = k END DO END DO ! On evite kupper = 1 et kupper = klev DO i = 1, klon kupper(i) = max(kupper(i), 2) kupper(i) = min(kupper(i), klev-1) END DO !---------- FIN nouveau calcul hw et ptop ------------------------------------- IF (iflag_wk_new_ptop==999) THEN DO i = 1, klon hw_(i)=0. ptop(i)=ph(i,1) Ktop(i)=1 pupper(i)=ph(i,2) kupper(i)=2 h_zzz(i)=0. Ptop1(i)=ph(i,1) ENDDO ENDIF zk_ptop_provis=k_ptop_provis RETURN END SUBROUTINE wake_pkupper END MODULE lmdz_wake_pkupper