Index: LMDZ6/trunk/libf/phylmd/coare_cp.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/coare_cp.F90	(revision 4722)
+++ LMDZ6/trunk/libf/phylmd/coare_cp.F90	(revision 5060)
@@ -1,52 +1,48 @@
-real function psit_30(zet)
-
-  real, intent(in) :: zet
-
-  if(zet<0) then
-     x=(1.-(15*zet))**.5 
-     psik=2*log((1+x)/2) 
-     x=(1.-(34.15*zet))**.3333 
-     psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 
-     f=zet*zet/(1+zet*zet) 
-     psit_30=(1-f)*psik+f*psic   
-
-  else
-     c=min(50.,.35*zet) 
-     psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
+module coare_cp_mod
+  implicit none
+  private
+  public psit_30, psiuo, coare_cp
+
+contains
+
+  real function psit_30(zet)
+    implicit none
+    real, intent(in) :: zet
+
+    real :: x, psik, psic, f, c
+
+    if(zet<0) then
+      x=(1.-(15.*zet))**.5
+      psik=2.*log((1.+x)/2)
+      x=(1.-(34.15*zet))**.3333
+      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
+      f=zet*zet/(1.+zet*zet)
+      psit_30=(1.-f)*psik+f*psic
+    else
+      c=min(50.,.35*zet)
+      psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
+    endif
+
+  end function psit_30
+
+  real function psiuo(zet)
+    implicit none
+    real, intent(in) :: zet
+
+    real :: x, psik, psic, f, c
+
+    if (zet<0) then
+       x=(1.-15.*zet)**.25
+       psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.)
+       x=(1.-10.15*zet)**.3333
+       psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
+       f=zet*zet/(1+zet*zet)
+       psiuo=(1-f)*psik+f*psic
+    else
+     c=min(50.,.35*zet)
+     psiuo=-((1.+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
   endif
-end FUNCTION psit_30
-
-real function psiuo(zet)
-
-  real, intent(in) :: zet
-
-  if (zet<0) then 
-
-     x=(1.-15.*zet)**.25 
-     psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) 
-     x=(1.-10.15*zet)**.3333 
-     psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 
-     f=zet*zet/(1+zet*zet) 
-     psiuo=(1-f)*psik+f*psic                                                
-  else
-     c=min(50.,.35*zet) 
-     psiuo=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
-  endif
-END FUNCTION psiuo
-
-
-real function gen_q1(ts,p,dq)
-
-  real, intent(in) :: ts,p,dq
-  real es,qs
-
-  es = 6.112*exp(17.502*(ts-273.15)/(ts-32.18))*.98*(1.0007+3.46e-8*p)
-  qs = es*62.197/(p-37.8*es)
-
-  q1 = dq + qs  
-
-end function gen_q1
-
-
+
+end function psiuo
 
 
@@ -61,5 +57,5 @@
   !version with shortened iteration  modified Rt and Rq
   !uses wave information wave period in s and wave ht in m
-  !no wave, standard coare 2.6 charnock:  jwave=0 
+  !no wave, standard coare 2.6 charnock:  jwave=0
   !Oost et al.  zo=50/2/pi L (u*/c)**4.5 if jwave=1
   !taylor and yelland  zo=1200 h*(L/h)**4.5 jwave=2
@@ -104,5 +100,5 @@
   real old_usr, old_tsr, old_qsr,tmp
 
-  real, external :: psit_30, psiuo, grv
+  real, external :: grv
   integer i,j,k
 !---------------- Rajout pour prendre en compte différent Z0 --------------------------------!
@@ -113,41 +109,41 @@
 !---------------------------------------------------------------------------------------------
 
-  Ribcu=-zu/zi/.004/Beta**3 
-  cpa=1004.67 
+  Ribcu=-zu/zi/.004/Beta**3
+  cpa=1004.67
 
   t_c = t - 273.15
 
 
-  Le=(2.501-.00237*(t_c-dt))*1e6 
-
-
-  visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) 
-
-
-  cpv=cpa*(1+0.84*Q) 
-  rhoa=P/(Rgas*t*(1+0.61*Q)) 
-
-
-
-
-
-
-  ug=.2 
-
-  ut=sqrt(du*du+ug*ug) 
+  Le=(2.501-.00237*(t_c-dt))*1e6
+
+
+  visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3)
+
+
+  cpv=cpa*(1+0.84*Q)
+  rhoa=P/(Rgas*t*(1+0.61*Q))
+
+
+
+
+
+
+  ug=.2
+
+  ut=sqrt(du*du+ug*ug)
   ut=MAX(ut , 0.1 * MIN(10.,zu) )
 
-  u10=ut*log(10/1e-4)/log(zu/1e-4) 
+  u10=ut*log(10/1e-4)/log(zu/1e-4)
   usr=.035*u10                              ! turbulent friction velocity (m/s), including gustiness
   zom10=0.011*usr*usr/grav+0.11*visa/usr     ! roughness length for u (smith 88)
-  Cd10=(von/log(10/zom10))**2 
+  Cd10=(von/log(10/zom10))**2
 !  zoh10=0.40*visa/usr+1.4e-5
-!  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca 
-  Ch10=0.00115 
-  Ct10=Ch10/sqrt(Cd10) 
+!  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca
+  Ch10=0.00115
+  Ct10=Ch10/sqrt(Cd10)
   zot10=10/exp(von/Ct10)                    ! roughness length for t
-  Cd=(von/log(zu/zom10))**2 
-  Ct=von/log(zt/zot10) 
-  CC=von*Ct/Cd 
+  Cd=(von/log(zu/zom10))**2
+  Ct=von/log(zt/zot10)
+  CC=von*Ct/Cd
 
   ut0 = ut
@@ -155,31 +151,31 @@
 
   ut = ut0
-  Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 
-
-  if (Ribu .LT. 0) then 
-     zetu=CC*Ribu/(1+Ribu/Ribcu) 
-  else 
-     zetu=CC*Ribu*(1+27/9*Ribu/CC) 
+  Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2
+
+  if (Ribu < 0) then
+     zetu=CC*Ribu/(1+Ribu/Ribcu)
+  else
+     zetu=CC*Ribu*(1+27/9*Ribu/CC)
   endif
 
-  L10=zu/zetu 
-
-  ! if (zetu .GT. 50) then 
-  !    nits=1 
+  L10=zu/zetu
+
+  ! if (zetu .GT. 50) then
+  !    nits=1
   ! endif
 
   usr=ut*von/(log(zu/zom10)-psiuo(zetu))
-  tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) 
-  qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) 
+  tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10))
+  qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10))
   !  tkt=.001
 
   ! charnock constant - lin par morceau - constant
   if ( ut <= 10. ) then
-     charn=0.011 
+     charn=0.011
   else
-     if (ut .GT. 18) then
+     if (ut > 18) then
         charn=0.018
      else
-        charn=0.011+(ut-10)/(18-10)*(0.018-0.011) 
+        charn=0.011+(ut-10)/(18-10)*(0.018-0.011)
      endif
   endif
@@ -192,5 +188,5 @@
 
   !***************  bulk loop ************
-  do i=1, nits 
+  do i=1, nits
 
      zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q)
@@ -200,17 +196,17 @@
 !     ELSE IF (NGRVWAVES==1) THEN
 !       zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 &
-!               + 0.11*visa/usr                       !Oost et al. 2002  
+!               + 0.11*visa/usr                       !Oost et al. 2002
 !     ELSE IF (NGRVWAVES==2) THEN
 !       zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 &
-!               + 0.11*visa/usr                       !Taulor and Yelland 2001  
-!     ENDIF 
-
-     zom=charn*usr*usr/grav+0.11*visa/usr 
-
-     rr=zom*usr/visa 
-     L=zu/zet 
+!               + 0.11*visa/usr                       !Taulor and Yelland 2001
+!     ENDIF
+
+     zom=charn*usr*usr/grav+0.11*visa/usr
+
+     rr=zom*usr/visa
+     L=zu/zet
      zoq=min(1.15e-4,5.5e-5/rr**.6)  ! a modifier
      zot=zoq                         ! a modifier
-!     zot=0.40*visa/usr+1.4e-5 
+!     zot=0.40*visa/usr+1.4e-5
 
      old_usr = usr
@@ -218,17 +214,17 @@
      old_qsr = qsr
 
-     usr=ut*von/(log(zu/zom)-psiuo(zu/L)) 
-     tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) 
-     qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) 
-
-     Bf=-grav/t*usr*(tsr+.61*t*qsr) 
-
-     if (Bf .GT. 0) then 
-        ug=Beta*(Bf*zi)**.333 
-     else 
-        ug=.2 
+     usr=ut*von/(log(zu/zom)-psiuo(zu/L))
+     tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L))
+     qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L))
+
+     Bf=-grav/t*usr*(tsr+.61*t*qsr)
+
+     if (Bf > 0) then
+        ug=Beta*(Bf*zi)**.333
+     else
+        ug=.2
      endif
 
-     ut=sqrt(du*du+ug*ug) 
+     ut=sqrt(du*du+ug*ug)
      ut=MAX(ut , 0.1 * MIN(10.,zu) )
 
@@ -237,6 +233,6 @@
 
   ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress
-  ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr 
-  ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr 
+  ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr
+  ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr
 
   tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du
@@ -249,4 +245,5 @@
   rugosh = zot
 
-
 end subroutine coare_cp
+
+end module coare_cp_mod
