Index: LMDZ5/branches/testing/libf/phylmd/cdrag.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cdrag.F90	(revision 2787)
+++ LMDZ5/branches/testing/libf/phylmd/cdrag.F90	(revision 2886)
@@ -2,22 +2,26 @@
 !$Id$
 !
- SUBROUTINE cdrag( knon,  nsrf,   &
+ SUBROUTINE cdrag(knon,  nsrf,   &
      speed, t1,    q1,    zgeop1, &
      psol,  tsurf, qsurf, z0m, z0h,  &
-     pcfm,  pcfh,  zri,   pref )
+     cdm,  cdh,  zri,   pref)
 
   USE dimphy
   USE indice_sol_mod
   USE print_control_mod, ONLY: lunout, prt_level
+  USE ioipsl_getin_p_mod, ONLY : getin_p
+
   IMPLICIT NONE
 ! ================================================================= c
 !
 ! Objet : calcul des cdrags pour le moment (pcfm) et 
-!         les flux de chaleur sensible et latente (pcfh).   
-!
-! Modified histroy: 
-!   27-Jan-2014: richardson number inconsistant between
-!   coefcdrag.F90 and clcdrag.F90, Fuxing WANG wrote this subroutine
-!   by merging coefcdrag and clcdrag.
+!         les flux de chaleur sensible et latente (pcfh) d'apr??s
+!         Louis 1982, Louis 1979, King et al 2001
+!         ou Zilitinkevich et al 2002  pour les cas stables, Louis 1979
+!         et 1982 pour les cas instables 
+!
+! Modified history: 
+!  writting on the 20/05/2016
+!  modified on the 13/12/2016 to be adapted to LMDZ6
 !
 ! References:
@@ -28,4 +32,7 @@
 !     parametrization, November 1981, ECMWF, Reading, England. 
 !     Page: 19. Equations in Table 1.
+!   Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS 
+!    IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM
+!    Boundary-Layer Meteorology 72: 331-344
 !   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer.  
 !     European Centre for Medium-Range Weather Forecasts.
@@ -34,103 +41,121 @@
 !     model to the parameterization of evaporation from the tropical oceans. J.
 !     Climate, 5:418-434.
+!   King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate
+!   to surface and boundary-layer flux parametrizations
+!   QJRMS, 127, pp 779-794
 !
 ! ================================================================= c
-!
-! knon----input-I- nombre de points pour un type de surface
-! nsrf----input-I- indice pour le type de surface; voir indicesol.h
-! speed---input-R- module du vent au 1er niveau du modele
-! t1------input-R- temperature de l'air au 1er niveau du modele
-! q1------input-R- humidite de l'air au 1er niveau du modele
-! zgeop---input-R- geopotentiel au 1er niveau du modele
-! psol----input-R- pression au sol
-! tsurf---input-R- temperature de l'air a la surface 
-! qsurf---input-R- humidite de l'air a la surface
-! z0m, z0h---input-R- rugosite
-!! u1, v1 are removed, speed is used. Fuxing WANG, 04/03/2015,
-!! u1------input-R- vent zonal au 1er niveau du modele
-!! v1------input-R- vent meridien au 1er niveau du modele
-!
-! pcfm---output-R- cdrag pour le moment 
-! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible
-! zri----output-R- Richardson number
-! pref---output-R- pression au niveau zgeop/RG
-!
-! Parameters: 
-! ckap-----Karman constant
-! cb,cc,cd-parameters in Louis et al., 1982
 ! ================================================================= c
-!
-!
+! On choisit le couple de fonctions de correction avec deux flags:
+! Un pour les cas instables, un autre pour les cas stables
+!
+! iflag_corr_insta:
+!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
+!                   2: Louis 1982 
+!                   3: Laurent Li
+!
+! iflag_corr_sta:
+!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
+!                   2: Louis 1982 
+!                   3: Laurent Li
+!                   4: King  2001 (SHARP)
+!                   5: MO 1st order theory (allow collapse of turbulence)
+!            
+!
+!*****************************************************************
 ! Parametres d'entree
 !***************************************************************** 
-  INTEGER, INTENT(IN)                      :: knon, nsrf
+  
+  INTEGER, INTENT(IN)                      :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
   REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
   REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
-  REAL, DIMENSION(klon), INTENT(IN)        :: psol  ! pression au sol
-  REAL, DIMENSION(klon), INTENT(IN)        :: t1    ! temperature at 1st level
-  REAL, DIMENSION(klon), INTENT(IN)        :: q1    ! humidity at 1st level
   REAL, DIMENSION(klon), INTENT(IN)        :: tsurf ! Surface temperature (K)
   REAL, DIMENSION(klon), INTENT(IN)        :: qsurf ! Surface humidity (Kg/Kg)
   REAL, DIMENSION(klon), INTENT(IN)        :: z0m, z0h ! Rugosity at surface (m) 
-!  paprs, pplay u1, v1: to be deleted
-!  they were in the old clcdrag. Fuxing WANG, 04/03/2015
-!  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
-!  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
-!  REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1
+  REAL, DIMENSION(klon), INTENT(IN)        :: t1  ! Temperature au premier niveau (K) 
+  REAL, DIMENSION(klon), INTENT(IN)        :: q1  ! humidite specifique au premier niveau (kg/kg)
+  REAL, DIMENSION(klon), INTENT(IN)        :: psol ! pression au sol
+
+
 
 ! Parametres de sortie
 !******************************************************************
-  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfm  ! Drag coefficient for heat flux
-  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfh  ! Drag coefficient for momentum
+  REAL, DIMENSION(klon), INTENT(OUT)       :: cdm  ! Drag coefficient for heat flux
+  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh  ! Drag coefficient for momentum
   REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
   REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
 
-! Parametres local
-  INTEGER                                  :: ng_q1    ! Number of grids that q1 < 0.0
-  INTEGER                                  :: ng_qsurf ! Number of grids that qsurf < 0.0
-!  zgeop1, psol: to be deleted, they are inputs now. Fuxing WANG, 04/03/2015
-!  REAL, DIMENSION(klon)                    :: zgeop1! geopotentiel au 1er niveau du modele
-!  REAL, DIMENSION(klon)                    :: psol  ! pression au sol
-!
-! ================================================================= c
-!
+! Variables Locales
+!******************************************************************
+ 
+
   INCLUDE "YOMCST.h"
   INCLUDE "YOETHF.h"
-!  INCLUDE "indicesol.h"
   INCLUDE "clesphys.h"
-!
-! Quelques constantes et options:
-!!$PB      REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
-  REAL, PARAMETER :: CKAP=0.40, CB=5.0, CC=5.0, CD=5.0, CEPDU2 = (0.1)**2
-!
-! Variables locales :
-  INTEGER               :: i
-  REAL                  :: zdu2, ztsolv
-  REAL                  :: ztvd, zscf
-  REAL                  :: zucf, zcr
-  REAL                  :: friv, frih
-  REAL, DIMENSION(klon) :: zcfm1, zcfm2 ! Drag coefficient for momentum
-  REAL, DIMENSION(klon) :: zcfh1, zcfh2 ! Drag coefficient for heat flux
-  LOGICAL, PARAMETER    :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li
-  REAL, DIMENSION(klon) :: zcdn_m, zcdn_h         ! Drag coefficient in neutral conditions
+
+
+  REAL, PARAMETER       :: CKAP=0.40, CKAPT=0.42
+  REAL                   CEPDU2
+  REAL                   ALPHA
+  REAL                   CB,CC,CD,C2,C3
+  REAL                   MU, CM, CH, B, CMstar, CHstar
+  REAL                   PM, PH, BPRIME
+  REAL                   C
+  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
+  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
+  INTEGER                i
+  REAL                   zdu2, ztsolv
+  REAL                   ztvd, zscf
+  REAL                   zucf, zcr
+  REAL                   friv, frih
+  REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
+  REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
   REAL zzzcd
-!
-! Fonctions thermodynamiques et fonctions d'instabilite
-  REAL                  :: fsta, fins, x
-  fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
-  fins(x) = SQRT(1.0-18.0*x)
+
+  LOGICAL, SAVE :: firstcall = .TRUE.
+  !$OMP THREADPRIVATE(firstcall)
+  INTEGER, SAVE :: iflag_corr_sta
+  !$OMP THREADPRIVATE(iflag_corr_sta)
+  INTEGER, SAVE :: iflag_corr_insta
+  !$OMP THREADPRIVATE(iflag_corr_insta)
+
+!===================================================================c
+! Valeurs numeriques des constantes
+!===================================================================c
+
+
+! Minimum du carre du vent
+
+ CEPDU2 = (0.1)**2
+
+! Louis 1982
+
+  CB=5.0
+  CC=5.0
+  CD=5.0
+
+
+! King 2001
+
+  C2=0.25
+  C3=0.0625
+
+  
+! Louis 1979
+
+  BPRIME=4.7
+  B=9.4
+  
+
+!MO
+
+  ALPHA=5.0
+  
 
 ! ================================================================= c
-! Fuxing WANG, 04/03/2015, delete the calculation of zgeop1 
-! (le geopotentiel du premier couche de modele).
-! zgeop1 is an input ivariable in this subroutine.
-!  DO i = 1, knon
-!     zgeop1(i) = RD * t1(i) / (0.5*(paprs(i,1)+pplay(i,1))) &
-!          * (paprs(i,1)-pplay(i,1))
-!  END DO
-! ================================================================= c
-!
+! Tests avant de commencer
 ! Fuxing WANG, 04/03/2015
 ! To check if there are negative q1, qsurf values.
+!====================================================================c
   ng_q1 = 0     ! Initialization
   ng_qsurf = 0  ! Initialization
@@ -154,80 +179,233 @@
   ENDIF
 
-! Calculer le frottement au sol (Cdrag)
-  DO i = 1, knon
-!------------------------------------------------------------
-! u1, v1 are replaced by speed. Fuxing WANG, 04/03/2015,
-!     zdu2 = MAX(CEPDU2,u1(i)**2+v1(i)**2)
-!------------------------------------------------------------
+
+
+!=============================================================================c
+! Calcul du cdrag
+!=============================================================================c
+
+! On choisit les fonctions de stabilite utilisees au premier appel
+!**************************************************************************
+  IF (firstcall) THEN
+   iflag_corr_sta=2
+   iflag_corr_insta=2
+ 
+   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
+   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
+
+   firstcall = .FALSE.
+ ENDIF
+
+!xxxxxxxxxxxxxxxxxxxxxxx
+  DO i = 1, knon        ! Boucle sur l'horizontal
+!xxxxxxxxxxxxxxxxxxxxxxx
+
+
+! calculs preliminaires:
+!***********************
+     
+
      zdu2 = MAX(CEPDU2, speed(i)**2)
-!     psol(i) = paprs(i,1)
      pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
-                 (1.+ RETV * max(q1(i),0.0))))  ! negative q1 set to zero
-!------------ the old calculations in clcdrag----------------
-!     ztsolv = tsurf(i) * (1.0+RETV*qsurf(i)) 
-!     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
-!          *(1.+RETV*q1(i)) 
-!------------------------------------------------------------
-! Fuxing WANG, 04/03/2015, in this revised version, 
-! the negative qsurf and q1 are set to zero (as in coefcdrag)
-     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero
-     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) &
-          *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero
+                 (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero 
+     ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0))    ! negative qsurf set to zero
+     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
+          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
      zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
 
 
-! Coefficients CD neutres pour m et h : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h))
+! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
+!********************************************************************
+
      zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
-     zcdn_m(i) = zzzcd*zzzcd
-     zcdn_h(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
-
-     IF (zri(i) .GT. 0.) THEN      ! situation stable
+     cdmn(i) = zzzcd*zzzcd
+     cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
+
+
+! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
+!*******************************************************
+
+!''''''''''''''
+! Cas instables
+!''''''''''''''
+
+ IF (zri(i) .LT. 0.) THEN    
+
+
+        SELECT CASE (iflag_corr_insta)
+    
+        CASE (1) ! Louis 1979 + Mascart 1995
+
+           MU=LOG(MAX(z0m(i)/z0h(i),0.01))
+           CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
+           PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
+           CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
+           PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
+           CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
+            & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
+            & * ((zgeop1(i)/(RG*z0h(i)))**PH)
+           CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
+            & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
+            & * ((zgeop1(i)/(RG*z0m(i)))**PM)
+         
+
+
+
+           FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
+           FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
+
+        CASE (2) ! Louis 1982
+
+           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
+                *(1.0+zgeop1(i)/(RG*z0m(i)))))
+           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
+           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
+
+
+        CASE (3) ! Laurent Li
+
+           
+           FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
+           FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
+
+
+
+         CASE default ! Louis 1982
+
+           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
+                *(1.0+zgeop1(i)/(RG*z0m(i)))))
+           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
+           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
+
+
+         END SELECT
+
+
+
+! Calcul des drags
+
+
+       cdm(i)=cdmn(i)*FM(i)
+       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
+
+
+! Traitement particulier des cas oceaniques
+! on applique Miller et al 1992 en l'absence de gustiness 
+
+  IF (nsrf == is_oce) THEN
+!        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
+
+        IF(iflag_gusts==0) THEN
+           zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
+           cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
+        ENDIF
+
+
+        cdm(i)=MIN(cdm(i),cdmmax)
+        cdh(i)=MIN(cdh(i),cdhmax)
+
+  END IF
+
+
+
+ ELSE                          
+
+!'''''''''''''''
+! Cas stables :
+!'''''''''''''''
         zri(i) = MIN(20.,zri(i))
-        IF (.NOT.zxli) THEN
+
+       SELECT CASE (iflag_corr_sta)
+    
+        CASE (1) ! Louis 1979 + Mascart 1995
+   
+           FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
+           FH(i)=FM(i)
+
+
+        CASE (2) ! Louis 1982
+           
            zscf = SQRT(1.+CD*ABS(zri(i)))
-           friv = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min)
-           zcfm1(i) = zcdn_m(i) * friv
-           frih = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min )
-!!$ PB     zcfh1(i) = zcdn(i) * frih
-!!$ PB     zcfh1(i) = f_cdrag_stable * zcdn(i) * frih
-           zcfh1(i) = f_cdrag_ter * zcdn_h(i) * frih
-           IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn_h(i) * frih
-!!$ PB
-           pcfm(i) = zcfm1(i)
-           pcfh(i) = zcfh1(i)
-        ELSE
-           pcfm(i) = zcdn_m(i)* fsta(zri(i))
-           pcfh(i) = zcdn_h(i)* fsta(zri(i))
-        ENDIF
-     ELSE                          ! situation instable
-        IF (.NOT.zxli) THEN
-           zucf = 1./(1.+3.0*CB*CC*zcdn_m(i)*SQRT(ABS(zri(i)) &
-                *(1.0+zgeop1(i)/(RG*z0m(i)))))
-           zcfm2(i) = zcdn_m(i)*amax1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
-!!$ PB     zcfh2(i) = zcdn_h(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
-           zcfh2(i) = f_cdrag_ter*zcdn_h(i)*amax1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
-           pcfm(i) = zcfm2(i)
-           pcfh(i) = zcfh2(i)
-        ELSE
-           pcfm(i) = zcdn_m(i)* fins(zri(i))
-           pcfh(i) = zcdn_h(i)* fins(zri(i))
-        ENDIF
-        IF(iflag_gusts==0) THEN
-! cdrah sur l'ocean cf. Miller et al. (1992) - only active when gustiness parameterization is not active
-           zcr = (0.0016/(zcdn_m(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
-           IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn_h(i)*(1.0+zcr**1.25)**(1./1.25)
-        ENDIF
-     ENDIF
-  END DO
-
+           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
+           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
+
+
+        CASE (3) ! Laurent Li
+   
+           FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
+           FH(i)=FM(i)
+
+
+        CASE (4)  ! King 2001
+          
+           if (zri(i) .LT. C2/2.) then
+           FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
+           FH(i)=  FM(i)
+
+
+           else
+           FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
+           FH(i)= FM(i)
+           endif
+
+
+        CASE (5) ! MO
+   
+          if (zri(i) .LT. 1./alpha) then
+
+           FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
+           FH(i)=FM(i)
+
+           else 
+
+
+           FM(i)=MAX(1E-7,f_ri_cd_min)
+           FH(i)=FM(i)
+
+          endif
+
+
+
+
+
+        CASE default ! Louis 1982
+
+           zscf = SQRT(1.+CD*ABS(zri(i)))
+           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
+           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
+
+
+
+   END SELECT
+
+! Calcul des drags
+
+
+       cdm(i)=cdmn(i)*FM(i)
+       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
+
+       IF(nsrf.EQ.is_oce) THEN
+
+        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
+        cdm(i)=MIN(cdm(i),cdmmax)
+        cdh(i)=MIN(cdh(i),cdhmax)
+
+      ENDIF
+
+
+ ENDIF
+
+
+
+
+!xxxxxxxxxxx
+  END DO   !  Fin de la boucle sur l'horizontal
+!xxxxxxxxxxx
 ! ================================================================= c
      
-  ! IM cf JLD : on seuille cdrag_m et cdrag_h
-  IF (nsrf == is_oce) THEN
-     DO i=1,knon
-        pcfm(i)=MIN(pcfm(i),cdmmax)
-        pcfh(i)=MIN(pcfh(i),cdhmax)
-     END DO
-  END IF
+ 
 
 END SUBROUTINE cdrag
+
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
