!
!$Id: cdrag_mod.f90 5285 2024-10-28 13:33:29Z abarral $
!
!
MODULE cdrag_mod
!
! This module contains some procedures for calculation of the cdrag
! coefficients for turbulent diffusion at surface
!
  IMPLICIT NONE

CONTAINS
!
!****************************************************************************************
!
!r original routine svn3623
!
 SUBROUTINE cdrag(knon,  nsrf,   &
     speed, t1,    q1,    zgeop1, &
     psol, pblh,  tsurf, qsurf, z0m, z0h,  &
     ri_in, iri_in, &
     cdm,  cdh,  zri,   pref, prain, tsol , pat1)

  USE dimphy
  USE coare_cp_mod, ONLY: coare_cp
  USE coare30_flux_cnrm_mod, ONLY: coare30_flux_cnrm
  USE indice_sol_mod
  USE print_control_mod, ONLY: lunout, prt_level
  USE ioipsl_getin_p_mod, ONLY : getin_p
  USE lmdz_atke_turbulence_ini, ONLY : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut
  USE yomcst_mod_h
  USE clesphys_mod_h
  USE yoethf_mod_h
  IMPLICIT NONE
! ================================================================= c
!
! Objet : calcul des cdrags pour le moment (pcfm) et 
!         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:
!   Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
!     atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202. 
!   Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
!     operational PBL parametrization at ECMWF'. Workshop on boundary layer
!     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.
!     Equations: 110-113. Page 40.
!   Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
!     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
! ================================================================= 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 ! 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)        :: tsurf    ! Surface temperature (K)
  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf    ! Surface humidity (Kg/Kg)
  REAL, DIMENSION(klon), INTENT(INOUT)     :: z0m, z0h ! Rugosity at surface (m) 
  REAL, DIMENSION(klon), INTENT(IN)        :: ri_in ! Input Richardson 1st layer for first guess calculations of screen var. 
  INTEGER, INTENT(IN)                      :: iri_in! iflag to activate cdrag calculation using ri1
  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

!------------------ Rajout pour COARE (OT2018) --------------------
  REAL, DIMENSION(klon), INTENT(INOUT)     :: pblh  !hauteur de CL
  REAL, DIMENSION(klon), INTENT(IN)        :: prain !rapport au precipitation
  REAL, DIMENSION(klon), INTENT(IN)        :: tsol  !SST imposé sur la surface oceanique
  REAL, DIMENSION(klon), INTENT(IN)        :: pat1  !pression premier lev



! Parametres de sortie
!******************************************************************
  REAL, DIMENSION(klon), INTENT(OUT)       :: cdm   ! Drag coefficient for momentum
  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh   ! Drag coefficient for heat flux
  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
  REAL, DIMENSION(klon), INTENT(INOUT)       :: pref  ! Pression au niveau zgeop/RG

! Variables Locales
!******************************************************************

  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
  INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
  INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
  INTEGER                i, k
  REAL                   zdu2, ztsolv
  REAL                   ztvd, zscf
  REAL                   zucf, zcr
  REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
  REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
  REAL zzzcd
  REAL, DIMENSION(klon) :: sm, prandtl              ! Stability function from atke turbulence scheme
  REAL                  ri0, ri1, cn                ! to have dimensionless term in sm and prandtl

!------------------ Rajout (OT2018) --------------------
!------------------ Rajout pour les appelles BULK (OT) --------------------
  REAL, DIMENSION(klon,2) :: rugos_itm 
  REAL, DIMENSION(klon,2) :: rugos_ith 
  REAL, PARAMETER         :: tol_it_rugos=1.e-4
  REAL, DIMENSION(3)      :: coeffs
  LOGICAL                 :: mixte
  REAL                    :: z_0m
  REAL                    :: z_0h
  REAL, DIMENSION(klon)   :: cdmm
  REAL, DIMENSION(klon)   :: cdhh

!------------------RAJOUT POUR ECUME -------------------

  REAL, DIMENSION(klon) :: PQSAT 
  REAL, DIMENSION(klon) :: PSFTH
  REAL, DIMENSION(klon) :: PFSTQ 
  REAL, DIMENSION(klon) :: PUSTAR
  REAL, DIMENSION(klon) :: PCD      ! Drag coefficient for momentum
  REAL, DIMENSION(klon) :: PCDN     ! Drag coefficient for momentum
  REAL, DIMENSION(klon) :: PCH      ! Drag coefficient for momentum
  REAL, DIMENSION(klon) :: PCE      ! Drag coefficient for momentum
  REAL, DIMENSION(klon) :: PRI 
  REAL, DIMENSION(klon) :: PRESA
  REAL, DIMENSION(klon) :: PSSS

  LOGICAL               :: OPRECIP
  LOGICAL               :: OPWEBB
  LOGICAL               :: OPERTFLUX
  LOGICAL               :: LPRECIP
  LOGICAL               :: LPWG



  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)
  LOGICAL, SAVE :: ok_cdrag_iter
  !$OMP THREADPRIVATE(ok_cdrag_iter)

!===================================================================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

! Consistent with atke scheme
  cn=(1./sqrt(cepsilon))**(2./3.)
  ri0=2./rpi*(cinf - cn)*ric/cn
  ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope


! ================================================================= 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
  DO i = 1, knon
     IF (q1(i).LT.0.0)     ng_q1 = ng_q1 + 1
     IF (qsurf(i).LT.0.0)  ng_qsurf = ng_qsurf + 1
  ENDDO
  IF (ng_q1.GT.0 .and. prt_level > 5) THEN
      WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"
      WRITE(lunout,*)" The total number of the grids is: ", ng_q1
      WRITE(lunout,*)" The negative q1 is set to zero "
!      abort_message="voir ci-dessus"
!      CALL abort_physic(modname,abort_message,1)
  ENDIF
  IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN
      WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"
      WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf
      WRITE(lunout,*)" The negative qsurf is set to zero "
!      abort_message="voir ci-dessus"
!      CALL abort_physic(modname,abort_message,1)
  ENDIF



!=============================================================================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
   ok_cdrag_iter = .FALSE.
 
   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
   CALL getin_p('ok_cdrag_iter',ok_cdrag_iter)

   firstcall = .FALSE.
 ENDIF

!------------------ Rajout (OT2018) --------------------
!---------    Rajout pour itération sur rugosité     ----------------
  rugos_itm(:,2) = z0m
  rugos_itm(:,1) = 3.*tol_it_rugos*z0m

  rugos_ith(:,2) = z0h !cp nouveau rugos_it
  rugos_ith(:,1) = 3.*tol_it_rugos*z0h
!--------------------------------------------------------------------

!xxxxxxxxxxxxxxxxxxxxxxx
  DO i = 1, knon        ! Boucle sur l'horizontal
!xxxxxxxxxxxxxxxxxxxxxxx


! calculs preliminaires:
!***********************
     
     zdu2 = MAX(CEPDU2, speed(i)**2)
     pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
          (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
     
!------------------ Rajout (OT2018) --------------------
!--------------   ON DUPLIQUE POUR METTRE ITERATION SUR OCEAN    ------------------------     
     IF (iri_in.EQ.1) THEN
       zri(i) = ri_in(i)
     ENDIF

     IF (nsrf == is_oce) THEN
        
!------------------  Pour Core 2 choix Core Pur ou Core Mixte  --------------------
       IF ((choix_bulk > 1 .AND. choix_bulk < 4) .AND. (nsrf .eq. is_oce)) THEN
         IF(choix_bulk .eq. 2) THEN
           mixte = .false.
         ELSE
            mixte = .true.
         ENDIF
         call clc_core_cp ( sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),t1(i),q1(i),&
             zgeop1(i)/RG, zgeop1(i)/RG,  zgeop1(i)/RG,&
             psol(i),nit_bulk,mixte,&
             coeffs,z_0m,z_0h)
         cdmm(i) = coeffs(1)
         cdhh(i) = coeffs(2)
         cdm(i)=cdmm(i)
         cdh(i)=cdhh(i)
         write(*,*) "clc_core cd ch",cdmm(i),cdhh(i)

!------------------                 Pour ECUME                 --------------------
       ELSE IF ((choix_bulk .eq. 4) .and. (nsrf .eq. is_oce)) THEN
         OPRECIP = .false.
         OPWEBB  = .false.
         OPERTFLUX = .false.
         IF (nsrf .eq. is_oce) THEN
           PSSS = 0.0
         ENDIF
         call ini_csts
         call ecumev6_flux( z_0m,t1(i),tsurf(i),&
             q1(i),qsurf(i),sqrt(zdu2),zgeop1(i)/RG,PSSS,zgeop1(i)/RG,&
             psol(i),pat1(i), OPRECIP, OPWEBB,&
             PSFTH,PFSTQ,PUSTAR,PCD,PCDN,&
             PCH,PCE,PRI,PRESA,prain,&
             z_0h,OPERTFLUX,coeffs)

         cdmm(i) = coeffs(1)
         cdhh(i) = coeffs(2)
         cdm(i)=cdmm(i)
         cdh(i)=cdhh(i)
   
         write(*,*) "ecume cd ch",cdmm(i),cdhh(i)

!------------------                 Pour COARE CNRM                 --------------------
       ELSE IF ((choix_bulk .eq. 5) .and. (nsrf .eq. is_oce)) THEN
         LPRECIP = .false.
         LPWG    = .false.
         call ini_csts
         block
           real, dimension(1) :: z0m_1d, z_0h_1d, sqrt_zdu2_1d, zgeop1_rg_1d  ! convert scalar to 1D for call
           z0m_1d = z0m
           z_0h_1d = z0h
           sqrt_zdu2_1d = sqrt(zdu2)
           zgeop1_rg_1d=zgeop1(i)/RG
           call coare30_flux_cnrm(z0m_1d,t1(i),tsurf(i), q1(i),  &
               sqrt_zdu2_1d,zgeop1_rg_1d,zgeop1_rg_1d,psol(i),qsurf(i),PQSAT, &
               PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, &
               PRESA,prain,pat1(i),z_0h_1d, LPRECIP, LPWG, coeffs)

         end block
         cdmm(i) = coeffs(1)
         cdhh(i) = coeffs(2)
         cdm(i)=cdmm(i)
         cdh(i)=cdhh(i)
         write(*,*) "coare CNRM cd ch",cdmm(i),cdhh(i)

!------------------                 Pour COARE Maison                 --------------------
       ELSE IF ((choix_bulk .eq. 1) .and. (nsrf .eq. is_oce)) THEN
         IF ( pblh(i) .eq. 0. ) THEN
           pblh(i) = 1500.
         ENDIF
         write(*,*) "debug size",size(coeffs)
         call coare_cp(sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),&
               t1(i),q1(i),&
               zgeop1(i)/RG,zgeop1(i)/RG,zgeop1(i)/RG,&
               psol(i), pblh(i),&
               nit_bulk,&
               coeffs,z_0m,z_0h)
         cdmm(i) = coeffs(1)
         cdhh(i) = coeffs(3)
         cdm(i)=cdmm(i)
         cdh(i)=cdhh(i)
         write(*,*) "coare cd ch",cdmm(i),cdhh(i)
       ELSE
!------------------                 Pour La param LMDZ (ocean)              --------------------
         zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
         IF (iri_in.EQ.1) THEN
           zri(i) = ri_in(i)
         ENDIF
       ENDIF
     

!----------------------- Rajout des itérations --------------
       DO  k=1,nit_bulk

! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
!********************************************************************
         zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*rugos_itm(i,2)))
         cdmn(i) = zzzcd*zzzcd
         cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*rugos_ith(i,2))))

! Calcul des fonctions de stabilite FMs, FHs, FMi, FHi :
!*******************************************************
         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 (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
                  sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
                  prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
                  FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
                  FH(i) = MAX((FM(i) / prandtl(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
           cdmm(i)=cdmn(i)*FM(i)
           cdhh(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.)
               cdhh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
             ENDIF
             cdmm(i)=MIN(cdmm(i),cdmmax)
             cdhh(i)=MIN(cdhh(i),cdhmax)
!             write(*,*) "LMDZ cd ch",cdmm(i),cdhh(i)
           END IF
         ELSE                          

!'''''''''''''''
! Cas stables :
!'''''''''''''''
           zri(i) = MIN(20.,zri(i))
           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)))
                  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 (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
                  sm(i) = MAX(smmin,cn*(1.-zri(i)/ric))
                           ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
                  prandtl(i) = pr_neut*exp(-pr_slope/pr_neut*zri(i)+zri(i)/pr_neut) &
                                + zri(i) * pr_slope
                  FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
             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
           
           cdmm(i)=cdmn(i)*FM(i)
           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
           
           IF (choix_bulk == 0) THEN
             cdm(i)=cdmn(i)*FM(i)
             cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
           ENDIF

           IF (nsrf.EQ.is_oce) THEN
             cdhh(i)=f_cdrag_oce*cdhn(i)*FH(i)
             cdmm(i)=MIN(cdmm(i),cdmmax)
             cdhh(i)=MIN(cdhh(i),cdhmax)
           ENDIF
           IF (ok_cdrag_iter) THEN
             rugos_itm(i,1) = rugos_itm(i,2)
             rugos_ith(i,1) = rugos_ith(i,2)
             rugos_itm(i,2) =  0.018*cdmm(i) * (speed(i))/RG  &
                              +  0.11*14e-6 / SQRT(cdmm(i) * zdu2)

!---------- Version SEPARATION DES Z0 ----------------------
             IF (iflag_z0_oce==0) THEN
               rugos_ith(i,2) = rugos_itm(i,2)
             ELSE IF (iflag_z0_oce==1) THEN
               rugos_ith(i,2) = 0.40*14e-6 / SQRT(cdmm(i) * zdu2)
             ENDIF
           ENDIF
         ENDIF
         IF (ok_cdrag_iter) THEN
           rugos_itm(i,2) = MAX(1.5e-05,rugos_itm(i,2))
           rugos_ith(i,2) = MAX(1.5e-05,rugos_ith(i,2))
         ENDIF
       ENDDO
       IF (nsrf.EQ.is_oce) THEN
         cdm(i)=MIN(cdmm(i),cdmmax)
         cdh(i)=MIN(cdhh(i),cdhmax)
       ENDIF
       z0m = rugos_itm(:,2)
       z0h = rugos_ith(:,2)
     ELSE ! (nsrf == is_oce)
       zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
       IF (iri_in.EQ.1) THEN
         zri(i) = ri_in(i)
       ENDIF

! 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)))
       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 (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
                 sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
                 prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
                 FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
                 FH(i) = MAX((FM(i) / prandtl(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)
       ELSE                          
!'''''''''''''''
! Cas stables :
!'''''''''''''''
         zri(i) = MIN(20.,zri(i))
         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)))
                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 (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
                sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
                prandtl(i) = pr_neut + zri(i) * pr_slope
                FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
                FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
           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)
       ENDIF
     ENDIF ! fin du if (nsrf == is_oce)
  END DO   !  Fin de la boucle sur l'horizontal

END SUBROUTINE cdrag

END MODULE cdrag_mod
