Ignore:
Timestamp:
Jul 16, 2024, 4:57:48 PM (2 months ago)
Author:
abarral
Message:

Refactor coare_cp.F90 & coare30_flux_cnrm.F90 into modules to reduce duplicated code. Accordingly adapt casting for coare30_flux_cnrm calls, which were inconsistent with declared interfaces.
Update deprecated math operators.
Lint end-of-line trailing spaces.
Remove dangling get_var1 from dynredem_mod.F90.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/coare_cp.F90

    r4722 r5060  
    1 real function psit_30(zet)
    2 
    3   real, intent(in) :: zet
    4 
    5   if(zet<0) then
    6      x=(1.-(15*zet))**.5
    7      psik=2*log((1+x)/2)
    8      x=(1.-(34.15*zet))**.3333
    9      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
    10      f=zet*zet/(1+zet*zet)
    11      psit_30=(1-f)*psik+f*psic   
    12 
    13   else
    14      c=min(50.,.35*zet)
    15      psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
     1module coare_cp_mod
     2  implicit none
     3  private
     4  public psit_30, psiuo, coare_cp
     5
     6contains
     7
     8  real function psit_30(zet)
     9    implicit none
     10    real, intent(in) :: zet
     11
     12    real :: x, psik, psic, f, c
     13
     14    if(zet<0) then
     15      x=(1.-(15.*zet))**.5
     16      psik=2.*log((1.+x)/2)
     17      x=(1.-(34.15*zet))**.3333
     18      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
     19      f=zet*zet/(1.+zet*zet)
     20      psit_30=(1.-f)*psik+f*psic
     21    else
     22      c=min(50.,.35*zet)
     23      psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
     24    endif
     25
     26  end function psit_30
     27
     28  real function psiuo(zet)
     29    implicit none
     30    real, intent(in) :: zet
     31
     32    real :: x, psik, psic, f, c
     33
     34    if (zet<0) then
     35       x=(1.-15.*zet)**.25
     36       psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.)
     37       x=(1.-10.15*zet)**.3333
     38       psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
     39       f=zet*zet/(1+zet*zet)
     40       psiuo=(1-f)*psik+f*psic
     41    else
     42     c=min(50.,.35*zet)
     43     psiuo=-((1.+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
    1644  endif
    17 end FUNCTION psit_30
    18 
    19 real function psiuo(zet)
    20 
    21   real, intent(in) :: zet
    22 
    23   if (zet<0) then
    24 
    25      x=(1.-15.*zet)**.25
    26      psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.)
    27      x=(1.-10.15*zet)**.3333
    28      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
    29      f=zet*zet/(1+zet*zet)
    30      psiuo=(1-f)*psik+f*psic                                               
    31   else
    32      c=min(50.,.35*zet)
    33      psiuo=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
    34   endif
    35 END FUNCTION psiuo
    36 
    37 
    38 real function gen_q1(ts,p,dq)
    39 
    40   real, intent(in) :: ts,p,dq
    41   real es,qs
    42 
    43   es = 6.112*exp(17.502*(ts-273.15)/(ts-32.18))*.98*(1.0007+3.46e-8*p)
    44   qs = es*62.197/(p-37.8*es)
    45 
    46   q1 = dq + qs 
    47 
    48 end function gen_q1
    49 
    50 
     45
     46end function psiuo
    5147
    5248
    … …  
    6157  !version with shortened iteration  modified Rt and Rq
    6258  !uses wave information wave period in s and wave ht in m
    63   !no wave, standard coare 2.6 charnock:  jwave=0 
     59  !no wave, standard coare 2.6 charnock:  jwave=0
    6460  !Oost et al.  zo=50/2/pi L (u*/c)**4.5 if jwave=1
    6561  !taylor and yelland  zo=1200 h*(L/h)**4.5 jwave=2
    … …  
    104100  real old_usr, old_tsr, old_qsr,tmp
    105101
    106   real, external :: psit_30, psiuo, grv
     102  real, external :: grv
    107103  integer i,j,k
    108104!---------------- Rajout pour prendre en compte différent Z0 --------------------------------!
    … …  
    113109!---------------------------------------------------------------------------------------------
    114110
    115   Ribcu=-zu/zi/.004/Beta**3 
    116   cpa=1004.67 
     111  Ribcu=-zu/zi/.004/Beta**3
     112  cpa=1004.67
    117113
    118114  t_c = t - 273.15
    119115
    120116
    121   Le=(2.501-.00237*(t_c-dt))*1e6 
    122 
    123 
    124   visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) 
    125 
    126 
    127   cpv=cpa*(1+0.84*Q) 
    128   rhoa=P/(Rgas*t*(1+0.61*Q)) 
    129 
    130 
    131 
    132 
    133 
    134 
    135   ug=.2 
    136 
    137   ut=sqrt(du*du+ug*ug) 
     117  Le=(2.501-.00237*(t_c-dt))*1e6
     118
     119
     120  visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3)
     121
     122
     123  cpv=cpa*(1+0.84*Q)
     124  rhoa=P/(Rgas*t*(1+0.61*Q))
     125
     126
     127
     128
     129
     130
     131  ug=.2
     132
     133  ut=sqrt(du*du+ug*ug)
    138134  ut=MAX(ut , 0.1 * MIN(10.,zu) )
    139135
    140   u10=ut*log(10/1e-4)/log(zu/1e-4) 
     136  u10=ut*log(10/1e-4)/log(zu/1e-4)
    141137  usr=.035*u10                              ! turbulent friction velocity (m/s), including gustiness
    142138  zom10=0.011*usr*usr/grav+0.11*visa/usr     ! roughness length for u (smith 88)
    143   Cd10=(von/log(10/zom10))**2 
     139  Cd10=(von/log(10/zom10))**2
    144140!  zoh10=0.40*visa/usr+1.4e-5
    145 !  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca 
    146   Ch10=0.00115 
    147   Ct10=Ch10/sqrt(Cd10) 
     141!  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca
     142  Ch10=0.00115
     143  Ct10=Ch10/sqrt(Cd10)
    148144  zot10=10/exp(von/Ct10)                    ! roughness length for t
    149   Cd=(von/log(zu/zom10))**2 
    150   Ct=von/log(zt/zot10) 
    151   CC=von*Ct/Cd 
     145  Cd=(von/log(zu/zom10))**2
     146  Ct=von/log(zt/zot10)
     147  CC=von*Ct/Cd
    152148
    153149  ut0 = ut
    … …  
    155151
    156152  ut = ut0
    157   Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 
    158 
    159   if (Ribu .LT. 0) then
    160      zetu=CC*Ribu/(1+Ribu/Ribcu) 
    161   else 
    162      zetu=CC*Ribu*(1+27/9*Ribu/CC) 
     153  Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2
     154
     155  if (Ribu < 0) then
     156     zetu=CC*Ribu/(1+Ribu/Ribcu)
     157  else
     158     zetu=CC*Ribu*(1+27/9*Ribu/CC)
    163159  endif
    164160
    165   L10=zu/zetu 
    166 
    167   ! if (zetu .GT. 50) then 
    168   !    nits=1 
     161  L10=zu/zetu
     162
     163  ! if (zetu .GT. 50) then
     164  !    nits=1
    169165  ! endif
    170166
    171167  usr=ut*von/(log(zu/zom10)-psiuo(zetu))
    172   tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) 
    173   qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) 
     168  tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10))
     169  qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10))
    174170  !  tkt=.001
    175171
    176172  ! charnock constant - lin par morceau - constant
    177173  if ( ut <= 10. ) then
    178      charn=0.011 
     174     charn=0.011
    179175  else
    180      if (ut .GT. 18) then
     176     if (ut > 18) then
    181177        charn=0.018
    182178     else
    183         charn=0.011+(ut-10)/(18-10)*(0.018-0.011) 
     179        charn=0.011+(ut-10)/(18-10)*(0.018-0.011)
    184180     endif
    185181  endif
    … …  
    192188
    193189  !***************  bulk loop ************
    194   do i=1, nits 
     190  do i=1, nits
    195191
    196192     zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q)
    … …  
    200196!     ELSE IF (NGRVWAVES==1) THEN
    201197!       zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 &
    202 !               + 0.11*visa/usr                       !Oost et al. 2002 
     198!               + 0.11*visa/usr                       !Oost et al. 2002
    203199!     ELSE IF (NGRVWAVES==2) THEN
    204200!       zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 &
    205 !               + 0.11*visa/usr                       !Taulor and Yelland 2001 
    206 !     ENDIF 
    207 
    208      zom=charn*usr*usr/grav+0.11*visa/usr 
    209 
    210      rr=zom*usr/visa 
    211      L=zu/zet 
     201!               + 0.11*visa/usr                       !Taulor and Yelland 2001
     202!     ENDIF
     203
     204     zom=charn*usr*usr/grav+0.11*visa/usr
     205
     206     rr=zom*usr/visa
     207     L=zu/zet
    212208     zoq=min(1.15e-4,5.5e-5/rr**.6)  ! a modifier
    213209     zot=zoq                         ! a modifier
    214 !     zot=0.40*visa/usr+1.4e-5 
     210!     zot=0.40*visa/usr+1.4e-5
    215211
    216212     old_usr = usr
    … …  
    218214     old_qsr = qsr
    219215
    220      usr=ut*von/(log(zu/zom)-psiuo(zu/L)) 
    221      tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) 
    222      qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) 
    223 
    224      Bf=-grav/t*usr*(tsr+.61*t*qsr) 
    225 
    226      if (Bf .GT. 0) then
    227         ug=Beta*(Bf*zi)**.333 
    228      else 
    229         ug=.2 
     216     usr=ut*von/(log(zu/zom)-psiuo(zu/L))
     217     tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L))
     218     qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L))
     219
     220     Bf=-grav/t*usr*(tsr+.61*t*qsr)
     221
     222     if (Bf > 0) then
     223        ug=Beta*(Bf*zi)**.333
     224     else
     225        ug=.2
    230226     endif
    231227
    232      ut=sqrt(du*du+ug*ug) 
     228     ut=sqrt(du*du+ug*ug)
    233229     ut=MAX(ut , 0.1 * MIN(10.,zu) )
    234230
    … …  
    237233
    238234  ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress
    239   ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr 
    240   ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr 
     235  ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr
     236  ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr
    241237
    242238  tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du
    … …  
    249245  rugosh = zot
    250246
    251 
    252247end subroutine coare_cp
     248
     249end module coare_cp_mod
Note: See TracChangeset for help on using the changeset viewer.