Ignore:
Timestamp:
Jan 16, 2020, 1:57:14 PM (5 years ago)
Author:
dubos
Message:

simple_physics : beautify code

Location:
dynamico_lmdz/simple_physics/phyparam
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • dynamico_lmdz/simple_physics/phyparam/Makefile

    r4228 r4229  
    2222clean :
    2323        rm -f obj/* include/* lib/* xml/*
    24 
     24        ../bash/beautify.sh physics/*.F90
    2525%.so : $(OBJECTS)
    2626        $(F90) -shared $^ -o $@

  • dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90

    r4228 r4229  
    11MODULE astronomy
    2  
     2
    33#include "use_logging.h"
    44
    … …  
    1212
    1313CONTAINS
    14  
     14
    1515  SUBROUTINE solarlong(pday,psollong)
    1616    REAL, INTENT(IN) :: pday               ! jour de l'annee (le jour 0 correspondant a l'equinoxe)
    … …  
    1818    LOGICAL, PARAMETER ::  lwrite=.TRUE.
    1919
    20 ! Local:
    21 ! ------
     20    ! Local:
     21    ! ------
    2222    REAL zanom,xref,zx0,zdx,zteta,zz
    2323    INTEGER iter
    2424
    25 !--------------------------------------------------------
    26 ! calcul de l'angle polaire et de la distance au soleil :
    27 ! -------------------------------------------------------
     25    !--------------------------------------------------------
     26    ! calcul de l'angle polaire et de la distance au soleil :
     27    ! -------------------------------------------------------
    2828
    29 !  calcul de l'zanomalie moyenne
     29    !  calcul de l'zanomalie moyenne
    3030
    3131    zz=(pday-peri_day)/year_day
    … …  
    3333    xref=abs(zanom)
    3434
    35 !  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
    36 !  methode de Newton
    37    
     35    !  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
     36    !  methode de Newton
     37
    3838    zx0=xref+e_elips*sin(xref)
    3939    DO iter=1,10
    … …  
    4343    zx0=zx0+zdx
    4444    if(zanom.lt.0.) zx0=-zx0
    45    
     45
    4646    ! zteta est la longitude solaire
    47    
     47
    4848    zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
    49    
     49
    5050    psollong=zteta-timeperi
    51    
     51
    5252    IF(psollong.LT.0.) psollong=psollong+2.*pi
    5353    IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
    … …  
    5555    !   sorties eventuelles:
    5656    !   ---------------------
    57    
     57
    5858    IF (lwrite) THEN
    5959       WRITELOG(*,*) 'day of the year  :',pday
    … …  
    6161       LOG_DBG('solarlong')
    6262    ENDIF
    63    
     63
    6464  END SUBROUTINE solarlong
    65  
     65
    6666  SUBROUTINE iniorbit
    6767    !=======================================================================
    … …  
    9292    !
    9393    !=======================================================================
    94    
     94
    9595    !-----------------------------------------------------------------------
    96    
     96
    9797    !   Local:
    9898    !   ------
    9999    REAL zxref,zanom,zz,zx0,zdx
    100100    INTEGER iter
    101    
     101
    102102    !-----------------------------------------------------------------------
    103    
     103
    104104    WRITELOG(*,*) 'Perihelie en Mkm  ',periheli
    105     WRITELOG(*,*) 'Aphelise  en Mkm  ',aphelie 
     105    WRITELOG(*,*) 'Aphelise  en Mkm  ',aphelie
    106106    WRITELOG(*,*) 'obliquite en degres  :',obliquit
    107    
     107
    108108    e_elips=(aphelie-periheli)/(periheli+aphelie)
    109109    p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
    110    
     110
    111111    WRITELOG(*,*) 'e_elips',e_elips
    112112    WRITELOG(*,*) 'p_elips',p_elips
    … …  
    115115    ! calcul de l'angle polaire et de la distance au soleil :
    116116    ! -------------------------------------------------------
    117    
     117
    118118    !  calcul de l'zanomalie moyenne
    119    
     119
    120120    zz=(year_day-peri_day)/year_day
    121121    zanom=2.*pi*(zz-nint(zz))
    122122    zxref=abs(zanom)
    123123    WRITELOG(*,*) 'zanom  ',zanom
    124    
     124
    125125    !  resolution de l'equation horaire  zx0 - e * sin (zx0) = zxref
    126126    !  methode de Newton
    127    
     127
    128128    zx0=zxref+e_elips*sin(zxref)
    129129    DO  iter=1,100
    130130       zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))
    131 !       if(abs(zdx).le.(1.e-12)) goto 120
    132131       zx0=zx0+zdx
    133132    END DO
    134    
     133
    135134    zx0=zx0+zdx
    136135    if(zanom.lt.0.) zx0=-zx0
    137136    WRITELOG(*,*) 'zx0   ',zx0
    138    
     137
    139138    ! zteta est la longitude solaire
    140    
     139
    141140    timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
    142141    WRITELOG(*,*) 'longitude solaire du perihelie timeperi = ',timeperi
    143    
     142
    144143    LOG_INFO('iniorbit')
    145    
     144
    146145  END SUBROUTINE iniorbit
    147146
    … …  
    171170    !   Declarations:
    172171    !   -------------
    173    
     172
    174173    ! arguments:
    175174    ! ----------
    … …  
    177176    REAL, INTENT(IN) :: pls
    178177    REAL, INTENT(OUT) :: pdist_sol,pdecli
    179    
     178
    180179    !-----------------------------------------------------------------------
    181    
     180
    182181    ! Distance Sun-Planet
    183    
     182
    184183    pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi))
    185    
     184
    186185    ! Solar declination
    187    
     186
    188187    pdecli= asin (sin(pls)*sin(obliquit*pi/180.))
    189    
     188
    190189  END SUBROUTINE orbite
    191190
    192191END MODULE astronomy
    193 

  • dynamico_lmdz/simple_physics/phyparam/physics/callkeys.F90

    r4223 r4229  
    22  IMPLICIT NONE
    33  SAVE
    4  
     4
    55  LOGICAL callrad,calldifv,calladj,callcond,callsoil,season,diurnal,lverbose
    66  INTEGER iradia

  • dynamico_lmdz/simple_physics/phyparam/physics/comgeomfi.F90

    r4217 r4229  
    55  REAL, ALLOCATABLE :: long(:), lati(:), sinlon(:), coslon(:), sinlat(:), coslat(:)
    66  INTEGER :: ngridmax, nlayermx, nsoilmx
    7 !$OMP THREADPRIVATE(long,lati,sinlon,coslon,sinlat,coslat,totarea)
    8 !$OMP THREADPRIVATE(ngridmax,nlayermx,nsoilmx)
     7  !$OMP THREADPRIVATE(long,lati,sinlon,coslon,sinlat,coslat,totarea)
     8  !$OMP THREADPRIVATE(ngridmax,nlayermx,nsoilmx)
    99
    1010CONTAINS
    11    
     11
    1212  SUBROUTINE init_comgeomfi(klon, klev, longitude, latitude)
    1313    INTEGER, INTENT(IN) :: klon, klev
    … …  
    2929    coslon(:)=cos(long(:))
    3030  END SUBROUTINE init_comgeomfi
    31  
     31
    3232END MODULE comgeomfi

  • dynamico_lmdz/simple_physics/phyparam/physics/convection.F90

    r4207 r4229  
    7676          ! * momentum u,v is mixed with weight zalpha, i.e. u:=zalpha*mean(u)+(1-zalpha)*u
    7777          ! where zalpha = sum( abs(theta-mean(theta))*dsig) / mean(theta)*sum(dsig)
    78           ! for large deviations of theta from its mean, this formula could in principe yield zalpha>1. 
     78          ! for large deviations of theta from its mean, this formula could in principe yield zalpha>1.
    7979          zalpha=0.
    8080          zum=0.
    … …  
    8484             zh2(i, l) = zhm
    8585             ! we must use compute zum, zvm from zu2, zv2 to conserve momentum (see below)
    86              zum=zum+dsig(l)*zu2(i,l) 
     86             zum=zum+dsig(l)*zu2(i,l)
    8787             zvm=zvm+dsig(l)*zv2(i,l)
    8888          END DO
    … …  
    101101             zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
    102102             zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
    103           ENDDO         
     103          ENDDO
    104104       END IF
    105        
     105
    106106    END DO
    107    
     107
    108108  END SUBROUTINE adjust_onecolumn
    109109
    … …  
    113113       &                  pdufi,pdvfi,pdhfi,     &
    114114       &                  pduadj,pdvadj,pdhadj)
    115    
     115
    116116    !=======================================================================
    117117    !
    … …  
    124124    !
    125125    !=======================================================================
    126    
     126
    127127    INTEGER, INTENT(IN) :: ngrid,nlay
    128128    REAL, INTENT(IN) ::  ptimestep
    … …  
    132132    REAL, INTENT(OUT) :: pdhadj(ngrid,nlay), pduadj(ngrid,nlay), pdvadj(ngrid,nlay)
    133133
    134     !   local:   
     134    !   local:
    135135    REAL sig(nlay+1),sdsig(nlay),dsig(nlay)
    136136    REAL zu(ngrid,nlay), zv(ngrid,nlay), zh(ngrid,nlay)
    137     REAL zu2(ngrid,nlay), zv2(ngrid,nlay), zh2(ngrid,nlay)   
     137    REAL zu2(ngrid,nlay), zv2(ngrid,nlay), zh2(ngrid,nlay)
    138138    LOGICAL vtest(ngrid)
    139139    INTEGER jadrs(ngrid)
    140140    INTEGER jcnt, ig, l, jj
    141    
     141
    142142    ! Add physics tendencies to u,v,h
    143143    DO l=1,nlay
    … …  
    151151    zv2(:,:)=zv(:,:)
    152152    zh2(:,:)=zh(:,:)
    153    
     153
    154154    !-----------------------------------------------------------------------
    155155    !   detection des profils a modifier:
    … …  
    160160    !   sinon, il reste a sa valeur initiale (.FALSE.)
    161161    !   cette operation est vectorisable
    162    
     162
    163163    vtest(:)=.FALSE.
    164164    DO l=2,nlay
    … …  
    186186       END DO
    187187    END DO
    188    
     188
    189189  END SUBROUTINE convadj
    190190

  • dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90

    r4228 r4229  
    99
    1010CONTAINS
    11  
     11
    1212  SUBROUTINE iniphyparam(ngrid,nlayer, &
    1313       &           punjours,              &
    … …  
    1818    USE planet, ONLY : coefir, coefvis
    1919    USE astronomy
    20     USE turbulence, ONLY : lmixmin, emin_turb 
     20    USE turbulence, ONLY : lmixmin, emin_turb
    2121    USE surface
    2222    USE read_param_mod
    … …  
    2727         pdayref        ! Day of reference for the simulation
    2828    REAL, INTENT(IN)  :: ptimestep, prad, pg, pr, pcpp, punjours
    29    
     29
    3030    CALL read_param('unjours', 86400., unjours,'unjours')
    3131    CALL read_param('planet_rad',prad,planet_rad,'planet_rad')
    … …  
    5353    CALL read_param('coefvis',.99     ,coefvis,'coefvis')
    5454    CALL read_param('coefir',.08      ,coefir,'coefir')
    55    
     55
    5656    CALL read_param('callrad',.true.,callrad,'appel rayonnemen')
    5757    CALL read_param('calldifv',.true.,calldifv,'appel difv')
    … …  
    6363    CALL read_param('lverbose',.true.,lverbose,'appel soil')
    6464    CALL read_param('period_sort',1.,period_sort,'period sorties en jour')
    65    
     65
    6666    CALL check_mismatch('unjours', punjours, unjours)
    6767    CALL check_mismatch('rad', prad, planet_rad)
    … …  
    6969    CALL check_mismatch('R', pr, r)
    7070    CALL check_mismatch('cpp', pcpp, cpp)
    71    
     71
    7272    WRITELOG(*,*) 'Activation de la physique:'
    7373    WRITELOG(*,*) ' R=',r

  • dynamico_lmdz/simple_physics/phyparam/physics/logging.F90

    r4227 r4229  
    1616     END SUBROUTINE plugin
    1717
    18      ! Plugin that writes into string 'line' information about the gridpoint of index 'index' 
     18     ! Plugin that writes into string 'line' information about the gridpoint of index 'index'
    1919     SUBROUTINE plugin_log_gridpoint(index, line)
    2020       INTEGER, INTENT(IN) :: index ! index of gridpoint
    … …  
    2525
    2626
    27   ! This module provides a default implementations of plugins but the top-level driver is welcome to override them. 
     27  ! This module provides a default implementations of plugins but the top-level driver is welcome to override them.
    2828  ! Note F2003/F2008: pgfortran (F2003) accepts to initialize pointers only to NULL()
    2929  ! => plugins are initialzed to NULL() and set to default values in flush_log and log_gridpoint
    … …  
    3939  INTEGER :: logging_lineno=0
    4040
    41   ! messages with a log level > max_log_level are not printed 
     41  ! messages with a log level > max_log_level are not printed
    4242  INTEGER, PARAMETER, PUBLIC :: log_level_fatal=1, log_level_error=2, log_level_warn=3, log_level_info=4, log_level_dbg=5
    4343  INTEGER, PUBLIC :: max_log_level = log_level_info
    … …  
    8181    logging_lineno = logging_lineno+1
    8282#ifndef XCODEML
    83     IF(.NOT.ASSOCIATED(log_gridpoint_plugin)) log_gridpoint_plugin => default_log_gridpoint 
     83    IF(.NOT.ASSOCIATED(log_gridpoint_plugin)) log_gridpoint_plugin => default_log_gridpoint
    8484    CALL log_gridpoint_plugin(index, logging_buf(logging_lineno))
    8585#endif
    … …  
    8787
    8888  SUBROUTINE default_log_gridpoint(index, line)
    89     INTEGER, INTENT(IN) :: index ! index of gridpoint                                                                                                                                   
     89    INTEGER, INTENT(IN) :: index ! index of gridpoint
    9090    CHARACTER(*), INTENT(OUT) :: line
    9191    line=''

  • dynamico_lmdz/simple_physics/phyparam/physics/mellor_yamada.F90

    r4205 r4229  
    55CONTAINS
    66
    7   PURE SUBROUTINE my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh) 
    8    
    9     !********************************************************************** 
    10     !****** FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) ******* 
    11     !* q2 au interfaces entre mailles.                                     
    12     !**********************************************************************
    13    
    14     INTEGER, INTENT(IN) :: imax,kmax 
     7  PURE SUBROUTINE my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh)
     8
     9    !**********************************************************************
     10    !****** FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) *******
     11    !* q2 au interfaces entre mailles.
     12    !**********************************************************************
     13
     14    INTEGER, INTENT(IN) :: imax,kmax
    1515    REAL, INTENT(IN) :: dt, gp
    1616    REAL, INTENT(IN) :: z(imax,kmax), &
    17          u(imax,kmax), v(imax,kmax), teta(imax,kmax), cd(imax)         
     17         u(imax,kmax), v(imax,kmax), teta(imax,kmax), cd(imax)
    1818    REAL, INTENT(INOUT) :: zi(imax,kmax+1), q2(imax,kmax+1)
    1919    REAL, INTENT(OUT) :: long(imax,kmax+1), km(imax,kmax+1), kh(imax,kmax+1)
    2020
    21     !************* DECLARATIONS ******************************************* 
    22    
    23     INTEGER, PARAMETER :: impmax=5 
     21    !************* DECLARATIONS *******************************************
     22
     23    INTEGER, PARAMETER :: impmax=5
    2424    REAL, PARAMETER :: kappa=0.4, &
    2525         a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08,           &
    26          e1=1.8, e2=1.33, &                                       
     26         e1=1.8, e2=1.33, &
    2727         khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5,              &
    28          q2min=0.001, q2lmin=0.001, &                             
    29          ghmax=0.023, ghmin=-0.28 
     28         q2min=0.001, q2lmin=0.001, &
     29         ghmax=0.023, ghmin=-0.28
    3030
    3131    REAL longc, au, ateta, az, adz, akq, acd, &
    … …  
    3333         am, azi, bet, beta, dest, du, dv, gh, &
    3434         prod, q2c, sh, sm, sq, us, vt, vt1, vt2
    35          
    36     REAL unsdz(imax,kmax),unsdzi(imax,kmax+1) 
     35
     36    REAL unsdz(imax,kmax),unsdzi(imax,kmax+1)
    3737    REAL kq(imax,kmax),                                               &
    38          &     m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1)             
     38         &     m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1)
    3939    REAL a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax),         &
    40          &     alph(imax,kmax)                                             
    41     REAL ksdz2inf,ksdz2sup 
     40         &     alph(imax,kmax)
     41    REAL ksdz2inf,ksdz2sup
    4242
    4343    INTEGER :: i,k
    44    
    45     !************* INCREMENTS VERTICAUX *********************************** 
    46    
    47     DO i=1,imax 
    48        zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax)) 
    49     END DO
    50     DO k=1,kmax 
    51        DO  i=1,imax 
    52           unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k)) 
    53        END DO
    54     END DO
    55 
    56     DO k=2,kmax 
    57        DO i=1,imax 
    58           unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1)) 
     44
     45    !************* INCREMENTS VERTICAUX ***********************************
     46
     47    DO i=1,imax
     48       zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax))
     49    END DO
     50    DO k=1,kmax
     51       DO  i=1,imax
     52          unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k))
     53       END DO
     54    END DO
     55
     56    DO k=2,kmax
     57       DO i=1,imax
     58          unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1))
    5959       END do
    6060    END DO
    6161
    62     DO i=1,imax 
    63        unsdzi(i,1)=0.5/(z(i,1)-zi(i,1)) 
    64        unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax)) 
     62    DO i=1,imax
     63       unsdzi(i,1)=0.5/(z(i,1)-zi(i,1))
     64       unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax))
    6565    END do
    66    
    67     !********************************************************************** 
    68                      
    69    
    70     !************* DIFFUSIVITES KH, KM et KQ ****************************** 
    71     ! Ci-dessous, une premiere estimation des diffusivites turbulentes km * 
    72     ! et kh est effectuee pour utilisation dans les taux de production    * 
    73     ! et de destruction de q2 et q2l. On calcule aussi kq.                * 
    74    
    75     DO k=2,kmax 
    76        DO i=1,imax 
    77           beta=2.0/(teta(i,k)+teta(i,k-1)) 
    78           n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1)) 
    79           n2(i,k)=amax1(0.0,n2(i,k)) 
    80           du=unsdzi(i,k)*(u(i,k)-u(i,k-1)) 
    81           dv=unsdzi(i,k)*(v(i,k)-v(i,k-1)) 
    82           m2(i,k)=du*du+dv*dv 
    83           ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10) 
    84           ri(i,k)=amax1(-0.1,min(4.0,ri(i,k))) 
    85        END DO
    86     END DO
    87 
    88     DO k=2,kmax 
    89        DO i=1,imax 
    90           vt=kappa*(zi(i,k)-zi(i,1)) 
    91           long(i,k)=vt/(1.0+vt/160.0) 
    92           IF(n2(i,k).gt.0.0) THEN 
    93              long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k)) 
     66
     67    !**********************************************************************
     68
     69
     70    !************* DIFFUSIVITES KH, KM et KQ ******************************
     71    ! Ci-dessous, une premiere estimation des diffusivites turbulentes km *
     72    ! et kh est effectuee pour utilisation dans les taux de production    *
     73    ! et de destruction de q2 et q2l. On calcule aussi kq.                *
     74
     75    DO k=2,kmax
     76       DO i=1,imax
     77          beta=2.0/(teta(i,k)+teta(i,k-1))
     78          n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1))
     79          n2(i,k)=amax1(0.0,n2(i,k))
     80          du=unsdzi(i,k)*(u(i,k)-u(i,k-1))
     81          dv=unsdzi(i,k)*(v(i,k)-v(i,k-1))
     82          m2(i,k)=du*du+dv*dv
     83          ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10)
     84          ri(i,k)=amax1(-0.1,min(4.0,ri(i,k)))
     85       END DO
     86    END DO
     87
     88    DO k=2,kmax
     89       DO i=1,imax
     90          vt=kappa*(zi(i,k)-zi(i,1))
     91          long(i,k)=vt/(1.0+vt/160.0)
     92          IF(n2(i,k).gt.0.0) THEN
     93             long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))
    9494          END IF
    9595          gh=amax1(ghmin,                                                 &
    96                &     min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))             
     96               &     min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
    9797          sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*                          &
    9898               &           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/     &
    99                &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))         
    100           km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) 
    101           sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) 
    102           kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) 
    103        END DO
    104     END DO
    105 
    106     DO i=1,imax 
    107        us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1))) 
    108        vt1=(b1**0.666667)*us*us 
     99               &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
     100          km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
     101          sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
     102          kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
     103       END DO
     104    END DO
     105
     106    DO i=1,imax
     107       us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1)))
     108       vt1=(b1**0.666667)*us*us
    109109       vt2=(b1**0.6666667)*kappa*kappa*                                 &
    110             &     m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1))                 
    111        !       q2(i,1)=amax1(q2min,vt1,vt2)                                   
    112        q2(i,1)=amax1(q2min,vt1) 
    113        long(i,1)=0.0 
    114        long(i,kmax+1)=long(i,kmax) 
    115        sq=0.2 
    116        kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq) 
    117     END DO
    118    
    119     DO k=2,kmax 
    120        DO i=1,imax 
    121           longc=0.5*(long(i,k)+long(i,k+1)) 
    122           q2c=0.5*(q2(i,k)+q2(i,k+1)) 
    123           sq=0.2 
    124           kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq) 
    125        END DO
    126     END DO
    127 
    128     !********************************************************************** 
    129                                                                        
    130    
    131     !************* CALCUL DE Q2 ******************************************* 
    132    
    133     DO k=2,kmax 
    134        DO i=1,imax 
    135           prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k))) 
     110            &     m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1))
     111       !       q2(i,1)=amax1(q2min,vt1,vt2)
     112       q2(i,1)=amax1(q2min,vt1)
     113       long(i,1)=0.0
     114       long(i,kmax+1)=long(i,kmax)
     115       sq=0.2
     116       kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq)
     117    END DO
     118
     119    DO k=2,kmax
     120       DO i=1,imax
     121          longc=0.5*(long(i,k)+long(i,k+1))
     122          q2c=0.5*(q2(i,k)+q2(i,k+1))
     123          sq=0.2
     124          kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq)
     125       END DO
     126    END DO
     127
     128    !**********************************************************************
     129
     130
     131    !************* CALCUL DE Q2 *******************************************
     132
     133    DO k=2,kmax
     134       DO i=1,imax
     135          prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k)))
    136136          dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k))                            &
    137                &           +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k)))                 
    138           IF(k.lt.kmax) THEN 
    139              ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k) 
    140           ELSE 
    141              ksdz2sup=0.0 
     137               &           +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k)))
     138          IF(k.lt.kmax) THEN
     139             ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k)
     140          ELSE
     141             ksdz2sup=0.0
    142142          END IF
    143           ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1) 
    144           b(i,k)=-ksdz2inf*dt 
    145           a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup) 
    146           c(i,k)=-ksdz2sup*dt 
    147           f(i,k)=q2(i,k)+dt*prod 
    148        END DO
    149     END DO
    150     DO i=1,imax 
     143          ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1)
     144          b(i,k)=-ksdz2inf*dt
     145          a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup)
     146          c(i,k)=-ksdz2sup*dt
     147          f(i,k)=q2(i,k)+dt*prod
     148       END DO
     149    END DO
     150    DO i=1,imax
    151151       f(i,2)=f(i,2)                                                    &
    152             &       +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1)                 
    153     END DO
    154    
    155     DO i=1,imax 
    156        alph(i,2)=a(i,2) 
    157     END DO
    158    
    159     DO k=3,kmax 
    160        DO i=1,imax 
    161           bet=b(i,k)/alph(i,k-1) 
    162           alph(i,k)=a(i,k)-bet*c(i,k-1) 
    163           f(i,k)=f(i,k)-bet*f(i,k-1) 
    164        END DO
    165     END DO
    166    
    167     DO i=1,imax 
    168        q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax)) 
    169        q2(i,kmax+1)=q2(i,kmax) 
    170     END DO
    171    
    172     DO k=kmax-1,2,-1 
    173        DO i=1,imax 
    174           q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k)) 
    175        END DO
    176     END DO
    177    
    178     DO i=1,imax 
    179        q2(i,2)=amax1(q2(i,2),q2(i,1)) 
    180     END DO
    181    
    182 !**********************************************************************
    183                                                                        
    184                                                                        
    185 !************* EVALUATION FINALE DE KH ET KM **************************
    186                                                                        
    187     DO k=2,kmax 
    188        DO i=1,imax 
    189           IF(n2(i,k).gt.0.0) THEN 
    190              long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))   
     152            &       +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1)
     153    END DO
     154
     155    DO i=1,imax
     156       alph(i,2)=a(i,2)
     157    END DO
     158
     159    DO k=3,kmax
     160       DO i=1,imax
     161          bet=b(i,k)/alph(i,k-1)
     162          alph(i,k)=a(i,k)-bet*c(i,k-1)
     163          f(i,k)=f(i,k)-bet*f(i,k-1)
     164       END DO
     165    END DO
     166
     167    DO i=1,imax
     168       q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax))
     169       q2(i,kmax+1)=q2(i,kmax)
     170    END DO
     171
     172    DO k=kmax-1,2,-1
     173       DO i=1,imax
     174          q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k))
     175       END DO
     176    END DO
     177
     178    DO i=1,imax
     179       q2(i,2)=amax1(q2(i,2),q2(i,1))
     180    END DO
     181
     182    !**********************************************************************
     183
     184
     185    !************* EVALUATION FINALE DE KH ET KM **************************
     186
     187    DO k=2,kmax
     188       DO i=1,imax
     189          IF(n2(i,k).gt.0.0) THEN
     190             long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))
    191191          END IF
    192192          gh=amax1(ghmin,                                                 &
    193                &         min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))         
     193               &         min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
    194194          sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*                          &
    195195               &           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/     &
    196                &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))         
    197           km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) 
    198           sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) 
    199           kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) 
    200        END DO
    201     END DO
    202    
    203     DO i=1,imax 
    204        km(i,1)=kmmin 
    205        km(i,kmax+1)=km(i,kmax) 
    206        kh(i,1)=khmin 
    207        kh(i,kmax+1)=kh(i,kmax) 
    208     END DO
    209    
    210     !********************************************************************** 
    211    
    212     am=1.0/float(imax) 
    213     DO k=kmax,1,-1 
    214        au=0.0 
    215        ateta=0.0 
    216        az=0.0 
    217        adz=0.0 
    218        akq=0.0 
    219        acd=0.0 
    220        DO i=1,imax 
    221           au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k)) 
    222           ateta=ateta+am*teta(i,k) 
    223           az=az+am*z(i,k) 
    224           adz=adz+am*(1.0/unsdz(i,k)) 
    225           akq=akq+am*kq(i,k) 
    226           acd=acd+am*cd(i) 
    227        END DO
    228     END DO
    229    
    230     DO k=kmax+1,1,-1 
    231        azi=0.0 
    232        adzi=0.0 
    233        aq2=0.0 
    234        al=0.0 
    235        akm=0.0 
    236        akh=0.0 
    237        am2=0.0 
    238        al0=0.0 
    239        DO i=1,imax 
    240           azi=azi+am*zi(i,k) 
    241           adzi=adzi+am*(1.0/unsdzi(i,k)) 
    242           aq2=aq2+am*q2(i,k) 
    243           al=al+am*long(i,k) 
    244           akm=akm+am*km(i,k) 
    245           akh=akh+am*kh(i,k) 
    246           am2=am2+am*m2(i,k) 
    247           !         al0=al0+am*long0d(i)                                         
    248        END DO
    249     END DO
    250    
     196               &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
     197          km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
     198          sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
     199          kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
     200       END DO
     201    END DO
     202
     203    DO i=1,imax
     204       km(i,1)=kmmin
     205       km(i,kmax+1)=km(i,kmax)
     206       kh(i,1)=khmin
     207       kh(i,kmax+1)=kh(i,kmax)
     208    END DO
     209
     210    !**********************************************************************
     211
     212    am=1.0/float(imax)
     213    DO k=kmax,1,-1
     214       au=0.0
     215       ateta=0.0
     216       az=0.0
     217       adz=0.0
     218       akq=0.0
     219       acd=0.0
     220       DO i=1,imax
     221          au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k))
     222          ateta=ateta+am*teta(i,k)
     223          az=az+am*z(i,k)
     224          adz=adz+am*(1.0/unsdz(i,k))
     225          akq=akq+am*kq(i,k)
     226          acd=acd+am*cd(i)
     227       END DO
     228    END DO
     229
     230    DO k=kmax+1,1,-1
     231       azi=0.0
     232       adzi=0.0
     233       aq2=0.0
     234       al=0.0
     235       akm=0.0
     236       akh=0.0
     237       am2=0.0
     238       al0=0.0
     239       DO i=1,imax
     240          azi=azi+am*zi(i,k)
     241          adzi=adzi+am*(1.0/unsdzi(i,k))
     242          aq2=aq2+am*q2(i,k)
     243          al=al+am*long(i,k)
     244          akm=akm+am*km(i,k)
     245          akh=akh+am*kh(i,k)
     246          am2=am2+am*m2(i,k)
     247          !         al0=al0+am*long0d(i)
     248       END DO
     249    END DO
     250
    251251  END SUBROUTINE my_25
    252  
     252
    253253END MODULE mellor_yamada

  • dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90

    r4228 r4229  
    1212       ps_rad=1.e5, height_scale=10000., ref_temp=285., &
    1313       capcal_nosoil=1e5, tsoil_init=300.
    14    
     14
    1515  ! internal state, written to / read from disk at checkpoint / restart
    1616  REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:), capcal(:), fluxgrd(:)
    … …  
    1919  ! precomputed arrays
    2020  REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
    21   !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)         
     21  !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)
    2222
    2323  INTEGER :: icount
    … …  
    4242
    4343    !=======================================================================
    44     !   Top routine of the physical parametrisations of the LMD 
     44    !   Top routine of the physical parametrisations of the LMD
    4545    !   20 parameters GCM for planetary atmospheres.
    4646    !   It includes:
    … …  
    5656         nlayer                   ! Number of vertical layers.
    5757    LOGICAL, INTENT(IN) ::      &
    58          firstcall,             & ! True at the first call 
     58         firstcall,             & ! True at the first call
    5959         lastcall                 ! True at the last call
    6060    REAL, INTENT(IN)    ::      &
    … …  
    7373         pdt(ngrid,nlayer),     &
    7474         pdpsrf(ngrid)
    75    
     75
    7676    !    Local variables :
    7777    REAL, DIMENSION(ngrid) :: mu0
    … …  
    9090    REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
    9191    REAL zfluxsw(ngrid),zfluxlw(ngrid), fluxrad(ngrid)
    92    
     92
    9393    WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid)
    9494    WRITELOG(*,*) 'nlayer',nlayer
    … …  
    104104
    105105    nlevel=nlayer+1
    106     igout=ngrid/2+1   
     106    igout=ngrid/2+1
    107107    zday=rjourvrai+gmtime
    108108
    … …  
    110110    !    0. Allocate and initialize at first call
    111111    !    --------------------
    112    
     112
    113113    IF(firstcall) THEN
    114114       !         zday_last=rjourvrai
    … …  
    122122    !    1. Initialisations :
    123123    !    --------------------
    124    
     124
    125125    icount=icount+1
    126    
     126
    127127    pdv(:,:)  = 0.
    128128    pdu(:,:)  = 0.
    … …  
    131131    zflubid(:)= 0.
    132132    zdtsrf(:) = 0.
    133    
     133
    134134    !-----------------------------------------------------------------------
    135135    !   calcul du geopotentiel aux niveaux intercouches
    136136    !   ponderation des altitudes au niveau des couches en dp/p
    137    
     137
    138138    DO l=1,nlayer
    139139       DO ig=1,ngrid
    … …  
    151151       ENDDO
    152152    ENDDO
    153    
     153
    154154    !-----------------------------------------------------------------------
    155155    !   Transformation de la temperature en temperature potentielle
    … …  
    160160       ENDDO
    161161    ENDDO
    162    
     162
    163163    !-----------------------------------------------------------------------
    164164    !    2. Compute radiative tendencies :
    165165    !    ---------------------------------------
    166    
     166
    167167    IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, &
    168             gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
    169             pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0)
     168         gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
     169         pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0)
    170170
    171171    !-----------------------------------------------------------------------
    … …  
    178178          zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
    179179       ENDDO
    180        
     180
    181181       zdum1(:,:)=0.
    182182       zdum2(:,:)=0.
    183        
     183
    184184       do l=1,nlayer
    185185          do ig=1,ngrid
    … …  
    187187          enddo
    188188       enddo
    189        
     189
    190190       CALL vdif(ngrid,nlayer,zday,        &
    191        &        ptimestep,capcal,z0,       &
    192        &        pplay,pplev,zzlay,zzlev,   &
    193        &        pu,pv,zh,tsurf,emissiv,    &
    194        &        zdum1,zdum2,zdum3,zflubid, &
    195        &        zdufr,zdvfr,zdhfr,zdtsrfr, &
    196        &        lverbose)
     191            &        ptimestep,capcal,z0,       &
     192            &        pplay,pplev,zzlay,zzlev,   &
     193            &        pu,pv,zh,tsurf,emissiv,    &
     194            &        zdum1,zdum2,zdum3,zflubid, &
     195            &        zdufr,zdvfr,zdhfr,zdtsrfr, &
     196            &        lverbose)
    197197
    198198       DO l=1,nlayer
    … …  
    203203          ENDDO
    204204       ENDDO
    205        
     205
    206206       DO ig=1,ngrid
    207207          zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
    208208       ENDDO
    209        
     209
    210210    ELSE
    211211       DO ig=1,ngrid
    212212          zdtsrf(ig)=zdtsrf(ig)+ &
    213           &      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
     213               &      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
    214214       ENDDO
    215215    ENDIF
    … …  
    218218    !   4. Dry convective adjustment:
    219219    !   -----------------------------
    220    
     220
    221221    IF(calladj) THEN
    222        
     222
    223223       DO l=1,nlayer
    224224          DO ig=1,ngrid
    … …  
    226226          ENDDO
    227227       ENDDO
    228        
     228
    229229       zdufr(:,:)=0.
    230230       zdvfr(:,:)=0.
    231231       zdhfr(:,:)=0.
    232        
     232
    233233       CALL convadj(ngrid,nlayer,ptimestep, &
    234        &                pplay,pplev,zpopsk, &
    235        &                pu,pv,zh,           &
    236        &                pdu,pdv,zdum1,      &
    237        &                zdufr,zdvfr,zdhfr)
    238        
     234            &                pplay,pplev,zpopsk, &
     235            &                pu,pv,zh,           &
     236            &                pdu,pdv,zdum1,      &
     237            &                zdufr,zdvfr,zdhfr)
     238
    239239       DO l=1,nlayer
    240240          DO ig=1,ngrid
    … …  
    244244          ENDDO
    245245       ENDDO
    246        
     246
    247247    ENDIF
    248248
    … …  
    250250    !   On ajoute les tendances physiques a la temperature du sol:
    251251    !   ---------------------------------------------------------------
    252    
     252
    253253    DO ig=1,ngrid
    254        tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) 
    255     ENDDO
    256    
     254       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
     255    ENDDO
     256
    257257    WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
    258    
     258
    259259    !-----------------------------------------------------------------------
    260260    !   soil temperatures:
    261261    !   --------------------
    262    
     262
    263263    IF (callsoil) THEN
    264264       CALL soil(ngrid,nsoilmx,.false.,inertie, &
    265        &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
     265            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
    266266       IF(lverbose) THEN
    267267          WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
    … …  
    271271       ENDIF
    272272    ENDIF
    273    
     273
    274274    !-----------------------------------------------------------------------
    275275    !   sorties:
    … …  
    285285       zday_last=zday
    286286       !  Ecriture/extension de la coordonnee temps
    287        
     287
    288288       do ig=1,ngridmax
    289289          zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*ref_temp))
    290290       enddo
    291        
     291
    292292       call writefield('u','Vent zonal moy','m/s',pu)
    293293       call writefield('v','Vent meridien moy','m/s',pv)
    … …  
    295295       call writefield('geop','Geopotential','m2/s2',pphi)
    296296       call writefield('plev','plev','Pa',pplev(:,1:nlayer))
    297        
     297
    298298       call writefield('du','du',' ',pdu)
    299299       call writefield('dv','du',' ',pdv)
    … …  
    301301       call writefield('dtsw','dtsw',' ',zdtsw)
    302302       call writefield('dtlw','dtlw',' ',zdtlw)
    303        
     303
    304304       call writefield('ts','Surface temper','K',tsurf)
    305305       call writefield('coslon','coslon',' ',coslon)
    … …  
    313313       call writefield('swsurf','SW surf','Pa',zfluxsw)
    314314       call writefield('lwsurf','LW surf','Pa',zfluxlw)
    315        
     315
    316316    endif
    317    
     317
    318318  END SUBROUTINE phyparam
    319319
    … …  
    338338         &                 zfluxsw(ngrid),      & ! short-wave flux at surface
    339339         &                 fluxrad(ngrid),      & ! surface flux
    340                            mu0(ngrid)             ! ??
     340         mu0(ngrid)             ! ??
    341341    ! local variables
    342342    REAL :: fract(ngrid), dtrad(ngrid, nlayer)
    … …  
    347347    !    2.1 Insolation
    348348    !    --------------------------------------------------
    349    
     349
    350350    CALL solarlong(zday,zls)
    351351    CALL orbite(zls,dist_sol,declin)
    352    
     352
    353353    IF(diurnal) THEN
    354354       IF ( .TRUE. ) then
    … …  
    362362          CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
    363363       ENDIF
    364        
     364
    365365       IF(lverbose) THEN
    366366          WRITELOG(*,*) 'day, declin, sinlon,coslon,sinlat,coslat'
    … …  
    375375       CALL mucorr(ngrid,declin,lati,mu0,fract,height_scale,planet_rad)
    376376    ENDIF
    377    
     377
    378378    zinsol=solarcst/(dist_sol*dist_sol)
    379    
     379
    380380    !    2.2 Radiative tendencies and fluxes:
    381381    !    --------------------------------------------------
    382    
     382
    383383    CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, &
    384384         &              pplev,ps_rad,                 &
    … …  
    386386         &              zfluxsw,zdtsw,                &
    387387         &              lverbose)
    388    
     388
    389389    CALL lw(ngrid,nlayer,coefir,emissiv, &
    390390         &             pplev,ps_rad,tsurf,pt, &
    391391         &             zfluxlw,zdtlw,         &
    392392         &             lverbose)
    393    
     393
    394394    !    2.4 surface fluxes
    395395    !    ------------------------------
    396    
     396
    397397    DO ig=1,ngrid
    398398       fluxrad(ig)=emissiv(ig)*zfluxlw(ig) &
    … …  
    402402       fluxrad(ig)=fluxrad(ig)-zplanck(ig)
    403403    ENDDO
    404    
     404
    405405    !    2.5 Temperature tendencies
    406406    !    --------------------------
    407    
     407
    408408    DO l=1,nlayer
    409409       DO ig=1,ngrid
    … …  
    412412       ENDDO
    413413    ENDDO
    414    
     414
    415415    IF(lverbose) THEN
    416416       WRITELOG(*,*) 'Diagnostics for radiation'
    … …  
    428428       LOG_DBG('radiative_tendencies')
    429429    ENDIF
    430    
     430
    431431  END SUBROUTINE radiative_tendencies
    432432

  • dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90

    r4215 r4229  
    11MODULE phys_const
    2  
     2
    33  REAL planet_rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg
    44

  • dynamico_lmdz/simple_physics/phyparam/physics/planet.F90

    r4190 r4229  
    44
    55  REAL :: coefvis, coefir
    6  
     6
    77END MODULE planet

  • dynamico_lmdz/simple_physics/phyparam/physics/radiative_lw.F90

    r4223 r4229  
    55  IMPLICIT NONE
    66  SAVE
    7  
     7
    88  PRIVATE
    9  
     9
    1010  PUBLIC :: lw
    11  
     11
    1212  LOGICAL, PARAMETER :: lstrong=.TRUE.
    1313  REAL, PARAMETER :: stephan=5.67e-08
    14  
     14
    1515CONTAINS
    16  
     16
    1717  SUBROUTINE lw(ngrid,nlayer,coefir,emissiv, &
    1818       pp,ps_rad,ptsurf,pt,              &
    … …  
    4545    !
    4646    !=======================================================================
    47    
     47
    4848    !   declarations:
    4949    !   -------------
    50    
     50
    5151    !   arguments:
    5252    !   ----------
    53    
     53
    5454    INTEGER, INTENT(IN)  ::  ngrid,nlayer
    5555    REAL,    INTENT(IN)  :: coefir,emissiv(ngrid),ps_rad
    … …  
    5757    REAL,    INTENT(OUT) ::  pdtlw(ngrid,nlayer),pfluxir(ngrid)
    5858    LOGICAL, INTENT(IN)  :: lwrite
    59    
     59
    6060    !   variables locales:
    6161    !   ------------------
    62    
     62
    6363    INTEGER nlevel,ilev,ig,i,il
    6464    REAL zplanck(ngrid,nlayer+1),zcoef
    … …  
    7272    !   initialisations:
    7373    !   ----------------
    74    
     74
    7575    nlevel=nlayer+1
    76    
     76
    7777    !-----------------------------------------------------------------------
    7878    !   2. calcul des quantites d'absorbants:
    7979    !   -------------------------------------
    80    
     80
    8181    !   absorption forte
    8282    IF(lstrong) THEN
    … …  
    9393       ENDIF
    9494       zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))
    95        
     95
    9696       !   absorption faible
    9797    ELSE
    … …  
    103103       zcoef=-log(coefir)/ps_rad
    104104    ENDIF
    105    
    106    
     105
     106
    107107    !-----------------------------------------------------------------------
    108108    !   2. calcul de la fonction de corps noir:
    109109    !   ---------------------------------------
    110    
     110
    111111    DO ilev=1,nlayer
    112112       DO ig=1,ngrid
    … …  
    116116       ENDDO
    117117    ENDDO
    118    
     118
    119119    !-----------------------------------------------------------------------
    120120    !   4. flux descendants:
    121121    !   --------------------
    122    
     122
    123123    DO ilev=1,nlayer
    124124       DO ig=1,ngrid
    … …  
    129129       ENDDO
    130130       CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
    131        
     131
    132132       DO il=nlayer,ilev,-1
    133133          zlwtr2(:)=zlwtr1(:)
    … …  
    142142       ENDDO
    143143    ENDDO
    144    
     144
    145145    DO ig=1,ngrid
    146146       zfluxdn(ig,nlevel)=0.
    147147       pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
    148148    ENDDO
    149    
     149
    150150    DO ig=1,ngrid
    151151       zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
    … …  
    153153            +(1.-emissiv(ig))*zfluxdn(ig,1)
    154154    ENDDO
    155    
     155
    156156    !-----------------------------------------------------------------------
    157157    !   3. flux montants:
    158158    !   ------------------
    159    
     159
    160160    DO ilev=1,nlayer
    161161       DO ig=1,ngrid
    … …  
    177177          ENDDO
    178178       ENDDO
    179        
     179
    180180    ENDDO
    181181
    … …  
    183183    !   5. calcul des flux nets:
    184184    !   ------------------------
    185    
     185
    186186    DO ilev=1,nlevel
    187187       DO ig=1,ngrid
    … …  
    189189       ENDDO
    190190    ENDDO
    191    
     191
    192192    !-----------------------------------------------------------------------
    193193    !   6. Calcul des taux de refroidissement:
    194194    !   --------------------------------------
    195    
     195
    196196    DO ilev=1,nlayer
    197197       DO ig=1,ngrid
    … …  
    200200       ENDDO
    201201    ENDDO
    202    
     202
    203203    !-----------------------------------------------------------------------
    204204    !   10. sorties eventuelles:
    205205    !   ------------------------
    206    
    207     IF (lwrite) THEN       
     206
     207    IF (lwrite) THEN
    208208       WRITELOG(*,*) 'Diagnostique rayonnement thermique'
    209209       WRITELOG(*,*) 'temperature     ', &
    210           'flux montant    flux desc.     taux de refroid.'
     210            'flux montant    flux desc.     taux de refroid.'
    211211       i=ngrid/2+1
    212212       WRITELOG(6,'(4e18.4)') ptsurf(i)
    … …  
    215215               zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
    216216       ENDDO
    217        WRITELOG(6,'(4e18.4)') zfluxup(i,nlevel),zfluxdn(i,nlevel)       
     217       WRITELOG(6,'(4e18.4)') zfluxup(i,nlevel),zfluxdn(i,nlevel)
    218218       LOG_DBG(tag)
    219219    ENDIF
    220220
    221 !-----------------------------------------------------------------------
    222    
     221    !-----------------------------------------------------------------------
     222
    223223  END SUBROUTINE lw
    224224
    … …  
    239239       ENDDO
    240240    ENDIF
    241    
     241
    242242  END SUBROUTINE lwtr
    243  
     243
    244244END MODULE radiative_lw

  • dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90

    r4228 r4229  
    55  IMPLICIT NONE
    66  SAVE
    7  
     7
    88  PRIVATE
    99
    … …  
    1111
    1212CONTAINS
    13  
    14   PURE SUBROUTINE monGATHER(n,a,b,index)     
     13
     14  PURE SUBROUTINE monGATHER(n,a,b,index)
    1515    INTEGER, INTENT(IN) ::  n,index(n)
    1616    REAL, INTENT(IN) :: b(n)
    … …  
    3232    END DO
    3333  end subroutine monscatter
    34  
     34
    3535  SUBROUTINE sw(ngrid,nlayer,ldiurn, &
    3636       coefvis,albedo, &
    … …  
    4141    !=======================================================================
    4242    !
    43     !   Rayonnement solaire en atmosphere non diffusante avec un 
     43    !   Rayonnement solaire en atmosphere non diffusante avec un
    4444    !   coefficient d'absoprption gris.
    4545    !
    … …  
    6060    REAL, INTENT(IN)    :: psolarf0
    6161    REAL, INTENT(OUT)   :: fsrfvis(ngrid),dtsw(ngrid,nlayer)
    62    
     62
    6363    REAL zalb(ngrid),zmu(ngrid),zfract(ngrid)
    6464    REAL zplev(ngrid,nlayer+1)
    6565    REAL zflux(ngrid),zdtsw(ngrid,nlayer)
    66    
     66
    6767    INTEGER ig,l,nlevel,index(ngrid),ncount,igout
    6868    REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1)
    … …  
    7171    REAL zu(ngrid,nlayer+1)
    7272    REAL tau0
    73    
     73
    7474    !-----------------------------------------------------------------------
    7575    !   1. initialisations:
    7676    !   -------------------
    77    
     77
    7878    nlevel=nlayer+1
    7979
    … …  
    8181    !   Definitions des tableaux locaux pour les points ensoleilles:
    8282    !   ------------------------------------------------------------
    83    
     83
    8484    IF (ldiurn) THEN
    8585       ncount=0
    … …  
    106106       zplev(:,:)=plevel(:,:)
    107107    ENDIF
    108    
     108
    109109    !-----------------------------------------------------------------------
    110110    !   calcul des profondeurs optiques integres depuis p=0:
    111111    !   ----------------------------------------------------
    112    
     112
    113113    tau0=-.5*log(coefvis)
    114    
     114
    115115    ! calcul de la partie homogene de l'opacite
    116116    tau0=tau0/ps_rad
    … …  
    120120       ENDDO
    121121    ENDDO
    122    
     122
    123123    !-----------------------------------------------------------------------
    124124    !   2. calcul de la transmission depuis le sommet de l'atmosphere:
    125125    !   -----------------------------------------------------------
    126    
     126
    127127    DO l=1,nlevel
    128128       DO ig=1,ncount
    … …  
    130130       ENDDO
    131131    ENDDO
    132    
     132
    133133    IF (lwrite) THEN
    134134       igout=ncount/2+1
    135        WRITELOG(*,*) 
     135       WRITELOG(*,*)
    136136       WRITELOG(*,*) 'Diagnostique des transmission dans le spectre solaire'
    137137       WRITELOG(*,*) 'zfract, zmu, zalb'
    … …  
    142142       ENDDO
    143143    ENDIF
    144    
     144
    145145    !-----------------------------------------------------------------------
    146146    !   3. taux de chauffage, ray. solaire direct:
    147147    !   ------------------------------------------
    148    
     148
    149149    DO l=1,nlayer
    150150       DO ig=1,ncount
    … …  
    155155    ENDDO
    156156    IF (lwrite) THEN
    157        WRITELOG(*,*) 
     157       WRITELOG(*,*)
    158158       WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
    159159       WRITELOG(*,*) ' 1 taux de chauffage lie au ray. solaire  direct'
    … …  
    162162       ENDDO
    163163    ENDIF
    164    
    165    
     164
     165
    166166    !-----------------------------------------------------------------------
    167167    !   4. calcul du flux solaire arrivant sur le sol:
    168168    !   ----------------------------------------------
    169    
     169
    170170    DO ig=1,ncount
    171171       z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1)
    … …  
    174174    ENDDO
    175175    IF (lwrite) THEN
    176        WRITELOG(*,*) 
     176       WRITELOG(*,*)
    177177       WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
    178178       WRITELOG(*,*) ' 2 flux solaire net incident sur le sol'
    179179       WRITELOG(*,*) zflux(igout)
    180180    ENDIF
    181    
    182    
     181
     182
    183183    !-----------------------------------------------------------------------
    184184    !   5.calcul des traansmissions depuis le sol, cas diffus:
    185185    !   ------------------------------------------------------
    186    
     186
    187187    DO l=1,nlevel
    188188       DO ig=1,ncount
    … …  
    190190       ENDDO
    191191    ENDDO
    192    
    193     IF (lwrite) THEN
    194        WRITELOG(*,*) 
     192
     193    IF (lwrite) THEN
     194       WRITELOG(*,*)
    195195       WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires'
    196196       WRITELOG(*,*) ' 3 transmission avec les sol'
    … …  
    200200       ENDDO
    201201    ENDIF
    202    
    203     !-----------------------------------------------------------------------
    204     !   6.ajout a l'echauffement de la contribution du ray. sol. reflechit: 
     202
     203    !-----------------------------------------------------------------------
     204    !   6.ajout a l'echauffement de la contribution du ray. sol. reflechit:
    205205    !   -------------------------------------------------------------------
    206    
     206
    207207    DO l=1,nlayer
    208208       DO ig=1,ncount
    … …  
    212212       ENDDO
    213213    ENDDO
    214    
     214
    215215    !-----------------------------------------------------------------------
    216216    !   10. sorites eventuelles:
    217217    !   ------------------------
    218    
    219     IF (lwrite) THEN
    220        WRITELOG(*,*) 
     218
     219    IF (lwrite) THEN
     220       WRITELOG(*,*)
    221221       WRITELOG(*,*) 'Diagnostique des taux de chauffage solaires:'
    222222       WRITELOG(*,*) ' 3 taux de chauffage total'
    … …  
    225225       ENDDO
    226226    ENDIF
    227    
     227
    228228    IF (ldiurn) THEN
    229229       fsrfvis(:)=0.
    … …  
    246246    !        call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir     ')
    247247    LOG_INFO('rad_sw')
    248    
     248
    249249  END SUBROUTINE sw
    250  
     250
    251251END MODULE radiative_sw
    252 

  • dynamico_lmdz/simple_physics/phyparam/physics/read_param_mod.F90

    r4227 r4229  
    55  SAVE
    66
    7   INTERFACE 
     7  INTERFACE
    88
    99     ! each of these plugins reads a parameter of a certain type and stores the value in 'val'
    … …  
    4343  INTERFACE read_param
    4444     PROCEDURE read_paramr, read_parami, read_paramb
    45   END INTERFACE read_param
     45  END INTERFACE
    4646
    4747  PUBLIC :: read_param
    … …  
    6868#endif
    6969  END SUBROUTINE read_parami
    70  
     70
    7171  SUBROUTINE read_paramb(name, defval, val, comment)
    7272    CHARACTER(*), INTENT(IN) :: name, comment
    … …  
    7878#endif
    7979  END SUBROUTINE read_paramb
    80  
     80
    8181END MODULE read_param_mod

  • dynamico_lmdz/simple_physics/phyparam/physics/solar.F90

    r4210 r4229  
    66  PRIVATE
    77
    8   REAL, PARAMETER :: pi_local = 4.0 * ATAN(1.0), deux_pi_local = 2.0 * pi_local 
     8  REAL, PARAMETER :: pi_local = 4.0 * ATAN(1.0), deux_pi_local = 2.0 * pi_local
    99
    1010  PUBLIC :: solang, zenang, mucorr
    … …  
    1414  PURE SUBROUTINE solang ( kgrid,psilon,pcolon,psilat,pcolat, &
    1515       ptim1,ptim2,ptim3,pmu0,pfract )
    16    
     16
    1717    !-----------------------------------------------------------------------
    1818    !          CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID
    … …  
    9393       ENDIF
    9494    ENDDO
    95    
     95
    9696  END SUBROUTINE solang
    9797
    9898  SUBROUTINE zenang(klon,longi,gmtime,pdtrad,lat,long,              &
    99        &                  pmu0,frac)                                     
    100     USE astronomy 
    101    
    102     !=============================================================         
    103     ! Auteur : O. Boucher (LMD/CNRS)                                       
    104     !          d'apres les routines zenith et angle de Z.X. Li             
    105     ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal     
    106     !          et l'ensoleillement moyen entre gmtime1 et gmtime2           
    107     !          connaissant la declinaison, la latitude et la longitude.     
    108     ! Rque   : Different de la routine angle en ce sens que zenang         
    109     !          fournit des moyennes de pmu0 et non des valeurs             
    110     !          instantanees, du coup frac prend toutes les valeurs         
    111     !          entre 0 et 1.                                               
    112     ! Date   : premiere version le 13 decembre 1994                         
    113     !          revu pour  GCM  le 30 septembre 1996                         
    114     !===============================================================       
    115     ! longi----INPUT : la longitude vraie de la terre dans son plan         
    116     !                  solaire a partir de l'equinoxe de printemps (degre) 
    117     ! gmtime---INPUT : temps universel en fraction de jour                 
    118     ! pdtrad---INPUT : pas de temps du rayonnement (secondes)               
    119     ! lat------INPUT : latitude en degres                                   
    120     ! long-----INPUT : longitude en degres                                 
    121     ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad   
    122     ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad   
    123     !================================================================       
    124     integer klon 
    125     !================================================================       
    126     real longi, gmtime, pdtrad 
    127     real lat(klon), long(klon), pmu0(klon), frac(klon) 
    128     !================================================================       
    129     integer i 
    130     real gmtime1, gmtime2 
    131     real incl 
    132     real omega1, omega2, omega 
    133     ! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi.       
    134     ! omega : heure en radian du coucher de soleil                         
    135     ! -omega est donc l'heure en radian de lever du soleil                 
    136     real omegadeb, omegafin 
    137     real zfrac1, zfrac2, z1_mu, z2_mu 
    138     ! declinaison en radian                     
    139     real lat_sun 
    140     ! longitude solaire en radian               
    141     real lon_sun 
    142     ! latitude du pt de grille en radian       
    143     real latr 
    144     !================================================================       
    145     !                                                                       
    146     !     incl=R_incl * pi_local / 180.                                     
     99       &                  pmu0,frac)
     100    USE astronomy
     101
     102    !=============================================================
     103    ! Auteur : O. Boucher (LMD/CNRS)
     104    !          d'apres les routines zenith et angle de Z.X. Li
     105    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
     106    !          et l'ensoleillement moyen entre gmtime1 et gmtime2
     107    !          connaissant la declinaison, la latitude et la longitude.
     108    ! Rque   : Different de la routine angle en ce sens que zenang
     109    !          fournit des moyennes de pmu0 et non des valeurs
     110    !          instantanees, du coup frac prend toutes les valeurs
     111    !          entre 0 et 1.
     112    ! Date   : premiere version le 13 decembre 1994
     113    !          revu pour  GCM  le 30 septembre 1996
     114    !===============================================================
     115    ! longi----INPUT : la longitude vraie de la terre dans son plan
     116    !                  solaire a partir de l'equinoxe de printemps (degre)
     117    ! gmtime---INPUT : temps universel en fraction de jour
     118    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
     119    ! lat------INPUT : latitude en degres
     120    ! long-----INPUT : longitude en degres
     121    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
     122    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
     123    !================================================================
     124    integer klon
     125    !================================================================
     126    real longi, gmtime, pdtrad
     127    real lat(klon), long(klon), pmu0(klon), frac(klon)
     128    !================================================================
     129    integer i
     130    real gmtime1, gmtime2
     131    real incl
     132    real omega1, omega2, omega
     133    ! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi.
     134    ! omega : heure en radian du coucher de soleil
     135    ! -omega est donc l'heure en radian de lever du soleil
     136    real omegadeb, omegafin
     137    real zfrac1, zfrac2, z1_mu, z2_mu
     138    ! declinaison en radian
     139    real lat_sun
     140    ! longitude solaire en radian
     141    real lon_sun
     142    ! latitude du pt de grille en radian
     143    real latr
     144    !================================================================
     145    !
     146    !     incl=R_incl * pi_local / 180.
    147147    WRITELOG(*,*) 'Obliquite =' ,obliquit
    148148    LOG_INFO('solar')
    149149
    150     incl=obliquit * pi_local / 180. 
    151     !                                                                       
    152     !     lon_sun = longi * pi_local / 180.0                               
    153     lon_sun = longi 
    154     lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) 
    155     !                                                                       
    156     gmtime1=gmtime*86400. 
    157     gmtime2=gmtime*86400.+pdtrad 
    158     !                                                                       
    159     DO i = 1, klon 
    160        !                                                                       
    161        !     latr = lat(i) * pi_local / 180.                                   
    162        latr = lat(i) 
    163        !                                                                       
    164        !--pose probleme quand lat=+/-90 degres                                 
    165        !                                                                       
    166        !      omega = -TAN(latr)*TAN(lat_sun)                                 
    167        !      omega = ACOS(omega)                                             
    168        !      IF (latr.GE.(pi_local/2.+lat_sun)                               
    169        !     .    .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN                   
    170        !         omega = 0.0       ! nuit polaire                             
    171        !      ENDIF                                                           
    172        !      IF (latr.GE.(pi_local/2.-lat_sun)                               
    173        !     .          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN             
    174        !         omega = pi_local  ! journee polaire                           
    175        !      ENDIF                                                           
    176        !                                                                       
    177        !--remplace par cela (le cas par defaut est different)                 
    178        !                                                                       
    179        !--nuit polaire                                       
    180        omega=0.0 
     150    incl=obliquit * pi_local / 180.
     151    !
     152    !     lon_sun = longi * pi_local / 180.0
     153    lon_sun = longi
     154    lat_sun = ASIN (SIN(lon_sun)*SIN(incl) )
     155    !
     156    gmtime1=gmtime*86400.
     157    gmtime2=gmtime*86400.+pdtrad
     158    !
     159    DO i = 1, klon
     160       !
     161       !     latr = lat(i) * pi_local / 180.
     162       latr = lat(i)
     163       !
     164       !--pose probleme quand lat=+/-90 degres
     165       !
     166       !      omega = -TAN(latr)*TAN(lat_sun)
     167       !      omega = ACOS(omega)
     168       !      IF (latr.GE.(pi_local/2.+lat_sun)
     169       !     .    .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN
     170       !         omega = 0.0       ! nuit polaire
     171       !      ENDIF
     172       !      IF (latr.GE.(pi_local/2.-lat_sun)
     173       !     .          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
     174       !         omega = pi_local  ! journee polaire
     175       !      ENDIF
     176       !
     177       !--remplace par cela (le cas par defaut est different)
     178       !
     179       !--nuit polaire
     180       omega=0.0
    181181       IF (latr.GE.(pi_local/2.-lat_sun)                                 &
    182             &          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN               
    183           ! journee polaire                           
    184           omega = pi_local 
     182            &          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
     183          ! journee polaire
     184          omega = pi_local
    185185       ENDIF
    186186       IF (latr.LT.(pi_local/2.+lat_sun).AND.                            &
    187187            &    latr.GT.(-pi_local/2.+lat_sun).AND.                           &
    188188            &    latr.LT.(pi_local/2.-lat_sun).AND.                            &
    189             &    latr.GT.(-pi_local/2.-lat_sun)) THEN                         
    190           omega = -TAN(latr)*TAN(lat_sun) 
    191           omega = ACOS(omega) 
    192        ENDIF
    193        !                                                                       
    194        omega1 = gmtime1 + long(i)*86400.0/360.0 
    195        omega1 = omega1 / 86400.0*deux_pi_local 
    196        omega1 = MOD (omega1+deux_pi_local, deux_pi_local) 
    197        omega1 = omega1 - pi_local 
    198        !                                                                       
    199        omega2 = gmtime2 + long(i)*86400.0/360.0 
    200        omega2 = omega2 / 86400.0*deux_pi_local 
    201        omega2 = MOD (omega2+deux_pi_local, deux_pi_local) 
    202        omega2 = omega2 - pi_local 
    203        !                                                                       
    204        !--on est dans la meme journee locale 
    205        IF (omega1.LE.omega2) THEN 
    206           !                                                                       
     189            &    latr.GT.(-pi_local/2.-lat_sun)) THEN
     190          omega = -TAN(latr)*TAN(lat_sun)
     191          omega = ACOS(omega)
     192       ENDIF
     193       !
     194       omega1 = gmtime1 + long(i)*86400.0/360.0
     195       omega1 = omega1 / 86400.0*deux_pi_local
     196       omega1 = MOD (omega1+deux_pi_local, deux_pi_local)
     197       omega1 = omega1 - pi_local
     198       !
     199       omega2 = gmtime2 + long(i)*86400.0/360.0
     200       omega2 = omega2 / 86400.0*deux_pi_local
     201       omega2 = MOD (omega2+deux_pi_local, deux_pi_local)
     202       omega2 = omega2 - pi_local
     203       !
     204       !--on est dans la meme journee locale
     205       IF (omega1.LE.omega2) THEN
     206          !
    207207          IF (omega2.LE.-omega .OR. omega1.GE.omega .OR. omega.LT.1e-5)  &
    208                THEN                                                           
    209              !--nuit           
    210              frac(i)=0.0 
    211              pmu0(i)=0.0 
     208               THEN
     209             !--nuit
     210             frac(i)=0.0
     211             pmu0(i)=0.0
    212212             !--jour+nuit/jou
    213           ELSE 
    214              omegadeb=MAX(-omega,omega1) 
    215              omegafin=MIN(omega,omega2) 
    216              frac(i)=(omegafin-omegadeb)/(omega2-omega1) 
     213          ELSE
     214             omegadeb=MAX(-omega,omega1)
     215             omegafin=MIN(omega,omega2)
     216             frac(i)=(omegafin-omegadeb)/(omega2-omega1)
    217217             pmu0(i)=SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*    &
    218                   (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)         
     218                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
    219219          ENDIF
    220           !                                                                       
    221           !---omega1 GT omega2 -- a cheval sur deux journees         
    222        ELSE 
    223           !                                                                       
    224           !-------------------entre omega1 et pi                                 
    225           !--nuit                               
    226           IF (omega1.GE.omega) THEN 
    227              zfrac1=0.0 
    228              z1_mu =0.0 
    229              !--jour+nuit                           
    230           ELSE 
    231              omegadeb=MAX(-omega,omega1) 
    232              omegafin=omega 
    233              zfrac1=omegafin-omegadeb 
     220          !
     221          !---omega1 GT omega2 -- a cheval sur deux journees
     222       ELSE
     223          !
     224          !-------------------entre omega1 et pi
     225          !--nuit
     226          IF (omega1.GE.omega) THEN
     227             zfrac1=0.0
     228             z1_mu =0.0
     229             !--jour+nuit
     230          ELSE
     231             omegadeb=MAX(-omega,omega1)
     232             omegafin=omega
     233             zfrac1=omegafin-omegadeb
    234234             z1_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*     &
    235                   (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)         
     235                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
    236236          ENDIF
    237           !---------------------entre -pi et omega2                               
    238           !--nuit                             
    239           IF (omega2.LE.-omega) THEN 
    240              zfrac2=0.0 
    241              z2_mu =0.0 
    242              !--jour+nuit                         
    243           ELSE 
    244              omegadeb=-omega 
    245              omegafin=MIN(omega,omega2) 
    246              zfrac2=omegafin-omegadeb 
     237          !---------------------entre -pi et omega2
     238          !--nuit
     239          IF (omega2.LE.-omega) THEN
     240             zfrac2=0.0
     241             z2_mu =0.0
     242             !--jour+nuit
     243          ELSE
     244             omegadeb=-omega
     245             omegafin=MIN(omega,omega2)
     246             zfrac2=omegafin-omegadeb
    247247             z2_mu =SIN(latr)*SIN(lat_sun) + COS(latr)*COS(lat_sun)*     &
    248                   (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)         
    249              !                                                                       
     248                  (SIN(omegafin)-SIN(omegadeb))/ (omegafin-omegadeb)
     249             !
    250250          ENDIF
    251           !-----------------------moyenne                                         
    252           frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1) 
    253           pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10) 
    254           !                                                                       
    255           !---comparaison omega1 et omega2                         
    256        ENDIF
    257        !                                                                       
     251          !-----------------------moyenne
     252          frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1)
     253          pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10)
     254          !
     255          !---comparaison omega1 et omega2
     256       ENDIF
     257       !
    258258    ENDDO
    259     !                                                                       
     259    !
    260260  END SUBROUTINE zenang
    261  
     261
    262262  PURE SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad)
    263    
     263
    264264    !=======================================================================
    265265    !
    266     !   Calcul of equivalent solar angle and and fraction of day whithout 
     266    !   Calcul of equivalent solar angle and and fraction of day whithout
    267267    !   diurnal cycle.
    268268    !
    … …  
    284284    !
    285285    !=======================================================================
    286    
     286
    287287    !-----------------------------------------------------------------------
    288288    !
    289289    !    0. Declarations :
    290290    !    -----------------
    291    
     291
    292292    !     Arguments :
    293293    !     -----------
    … …  
    306306
    307307    !-----------------------------------------------------------------------
    308    
     308
    309309    z = pdeclin
    310310    cz = cos (z)
    311311    sz = sin (z)
    312    
     312
    313313    DO j = 1, npts
    314        
     314
    315315       phi = plat(j)
    316316       cphi = cos(phi)
    … …  
    322322       a = 1.0 - t*t
    323323       ap = a
    324        
     324
    325325       IF(t.eq.0.) then
    326326          t=0.5*pi
    … …  
    334334          ENDIF
    335335       ENDIF
    336        
     336
    337337       pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi
    338338       pfract(j) = t / pi
    … …  
    344344       pmu(j) = pmu(j) / pfract(j)
    345345       IF (pmu(j).eq.0.) pfract(j) = 0.
    346        
     346
    347347    END DO
    348    
     348
    349349    !-----------------------------------------------------------------------
    350350    !   correction de rotondite:
    351351    !   ------------------------
    352    
     352
    353353    alph=phaut/prad
    354354    DO j=1,npts
    … …  
    357357       !    $          (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j))
    358358    END DO
    359    
     359
    360360  END SUBROUTINE mucorr
    361361

  • dynamico_lmdz/simple_physics/phyparam/physics/surface.F90

    r4228 r4229  
    1414
    1515  ! precomputed variables
    16   REAL :: lambda 
     16  REAL :: lambda
    1717  REAL,ALLOCATABLE :: dz1(:),dz2(:)
    1818  !$OMP THREADPRIVATE(dz1,dz2,zc,lambda)
    1919
    2020  ! internal state variables needed for checkpoint/restart
    21   REAL, ALLOCATABLE :: zc(:,:),zd(:,:) 
     21  REAL, ALLOCATABLE :: zc(:,:),zd(:,:)
    2222  !$OMP THREADPRIVATE(zc,zd)
    2323
    … …  
    2828  SUBROUTINE init_soil(ngrid,nsoil)
    2929    INTEGER, INTENT(IN) :: ngrid, nsoil
    30     REAL min_period,dalph_soil 
    31     REAL fz,rk,fz1,rk1,rk2 
     30    REAL min_period,dalph_soil
     31    REAL fz,rk,fz1,rk1,rk2
    3232    INTEGER :: jk
    3333
    3434    ! this is a function definition
    35     fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) 
    36    
    37     !-----------------------------------------------------------------------
    38     !   ground levels                                                       
    39     !   grnd=z/l where l is the skin depth of the diurnal cycle:           
    40     !   --------------------------------------------------------           
    41    
     35    fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
     36
     37    !-----------------------------------------------------------------------
     38    !   ground levels
     39    !   grnd=z/l where l is the skin depth of the diurnal cycle:
     40    !   --------------------------------------------------------
     41
    4242    WRITELOG(*,*) 'nsoil,ngrid,firstcall=',nsoil,ngrid, .TRUE.
    4343
    44     ALLOCATE(dz1(nsoil),dz2(nsoil)) 
    45    
    46     min_period=20000. 
    47     dalph_soil=2. 
    48    
    49     !   la premiere couche represente un dixieme de cycle diurne           
    50     fz1=sqrt(min_period/pi) 
    51    
    52     DO jk=1,nsoil 
    53        rk1=jk 
    54        rk2=jk-1 
    55        dz2(jk)=fz(rk1)-fz(rk2) 
    56     ENDDO
    57     DO jk=1,nsoil-1 
    58        rk1=jk+.5 
    59        rk2=jk-.5 
    60        dz1(jk)=1./(fz(rk1)-fz(rk2)) 
    61     ENDDO
    62     lambda=fz(.5)*dz1(1) 
    63 
    64     WRITELOG(*,*) 'full layers, intermediate layers (secoonds)' 
    65     DO jk=1,nsoil 
    66        rk=jk 
    67        rk1=jk+.5 
    68        rk2=jk-.5 
     44    ALLOCATE(dz1(nsoil),dz2(nsoil))
     45
     46    min_period=20000.
     47    dalph_soil=2.
     48
     49    !   la premiere couche represente un dixieme de cycle diurne
     50    fz1=sqrt(min_period/pi)
     51
     52    DO jk=1,nsoil
     53       rk1=jk
     54       rk2=jk-1
     55       dz2(jk)=fz(rk1)-fz(rk2)
     56    ENDDO
     57    DO jk=1,nsoil-1
     58       rk1=jk+.5
     59       rk2=jk-.5
     60       dz1(jk)=1./(fz(rk1)-fz(rk2))
     61    ENDDO
     62    lambda=fz(.5)*dz1(1)
     63
     64    WRITELOG(*,*) 'full layers, intermediate layers (secoonds)'
     65    DO jk=1,nsoil
     66       rk=jk
     67       rk1=jk+.5
     68       rk2=jk-.5
    6969       WRITELOG(*,*) fz(rk1)*fz(rk2)*pi,        &
    70             &        fz(rk)*fz(rk)*pi                                         
     70            &        fz(rk)*fz(rk)*pi
    7171    ENDDO
    7272    LOG_INFO('init_soil')
    7373
    7474  END SUBROUTINE init_soil
    75  
     75
    7676  SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i,          &
    7777       &          ptimestep,ptsrf,ptsoil,                  &
    78        &          pcapcal,pfluxgrd)                                       
    79    
     78       &          pcapcal,pfluxgrd)
     79
    8080    !=======================================================================
    81     !                                                                       
    82     !   Auteur:  Frederic Hourdin     30/01/92                             
    83     !   -------                                                             
    84     !                                                                       
    85     !   objet:  computation of : the soil temperature evolution             
    86     !   ------                   the surfacic heat capacity "Capcal"       
    87     !                            the surface conduction flux pcapcal       
    88     !                                                                       
    89     !                                                                       
    90     !   Method: implicit time integration                                   
    91     !   -------                                                             
    92     !   Consecutive ground temperatures are related by:                     
    93     !           T(k+1) = C(k) + D(k)*T(k)  (1)                             
    94     !   the coefficients C and D are computed at the t-dt time-step.       
    95     !   Routine structure:                                                 
    96     !   1)new temperatures are computed  using (1)                         
    97     !   2)C and D coefficients are computed from the new temperature       
    98     !     profile for the t+dt time-step                                   
    99     !   3)the coefficients A and B are computed where the diffusive         
    100     !     fluxes at the t+dt time-step is given by                         
    101     !            Fdiff = A + B Ts(t+dt)                                     
    102     !     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt                   
    103     !            with F0 = A + B (Ts(t))                                   
    104     !                 Capcal = B*dt                                         
    105     !                                                                       
    106     !   Interface:                                                         
    107     !   ----------                                                         
    108     !                                                                       
    109     !   Arguments:                                                         
    110     !   ----------                                                         
    111     !   ngrid               number of grid-points                           
    112     !   ptimestep              physical timestep (s)                       
    113     !   pto(ngrid,nsoil)     temperature at time-step t (K)                 
    114     !   ptn(ngrid,nsoil)     temperature at time step t+dt (K)             
    115     !   pcapcal(ngrid)      specific heat (W*m-2*s*K-1)                     
    116     !   pfluxgrd(ngrid)      surface diffusive flux from ground (Wm-2)     
    117     !                                                                       
     81    !
     82    !   Auteur:  Frederic Hourdin     30/01/92
     83    !   -------
     84    !
     85    !   objet:  computation of : the soil temperature evolution
     86    !   ------                   the surfacic heat capacity "Capcal"
     87    !                            the surface conduction flux pcapcal
     88    !
     89    !
     90    !   Method: implicit time integration
     91    !   -------
     92    !   Consecutive ground temperatures are related by:
     93    !           T(k+1) = C(k) + D(k)*T(k)  (1)
     94    !   the coefficients C and D are computed at the t-dt time-step.
     95    !   Routine structure:
     96    !   1)new temperatures are computed  using (1)
     97    !   2)C and D coefficients are computed from the new temperature
     98    !     profile for the t+dt time-step
     99    !   3)the coefficients A and B are computed where the diffusive
     100    !     fluxes at the t+dt time-step is given by
     101    !            Fdiff = A + B Ts(t+dt)
     102    !     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
     103    !            with F0 = A + B (Ts(t))
     104    !                 Capcal = B*dt
     105    !
     106    !   Interface:
     107    !   ----------
     108    !
     109    !   Arguments:
     110    !   ----------
     111    !   ngrid               number of grid-points
     112    !   ptimestep              physical timestep (s)
     113    !   pto(ngrid,nsoil)     temperature at time-step t (K)
     114    !   ptn(ngrid,nsoil)     temperature at time step t+dt (K)
     115    !   pcapcal(ngrid)      specific heat (W*m-2*s*K-1)
     116    !   pfluxgrd(ngrid)      surface diffusive flux from ground (Wm-2)
     117    !
    118118    !=======================================================================
    119     !   declarations:                                                       
    120     !   -------------                                                       
    121    
    122    
    123     !-----------------------------------------------------------------------
    124     !  arguments                                                           
    125     !  ---------                                                           
    126    
    127     INTEGER ngrid,nsoil 
    128     REAL ptimestep 
    129     REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid) 
    130     REAL pcapcal(ngrid),pfluxgrd(ngrid) 
    131     LOGICAL firstcall 
    132    
    133    
    134     !-----------------------------------------------------------------------
    135     !  local arrays                                                         
    136     !  ------------                                                         
    137    
    138     INTEGER ig,jk 
    139     REAL zdz2(nsoil),z1(ngrid) 
    140        
    141     IF (firstcall) THEN 
     119    !   declarations:
     120    !   -------------
     121
     122
     123    !-----------------------------------------------------------------------
     124    !  arguments
     125    !  ---------
     126
     127    INTEGER ngrid,nsoil
     128    REAL ptimestep
     129    REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid)
     130    REAL pcapcal(ngrid),pfluxgrd(ngrid)
     131    LOGICAL firstcall
     132
     133
     134    !-----------------------------------------------------------------------
     135    !  local arrays
     136    !  ------------
     137
     138    INTEGER ig,jk
     139    REAL zdz2(nsoil),z1(ngrid)
     140
     141    IF (firstcall) THEN
    142142       CALL init_soil(ngrid, nsoil)
    143143    ELSE
    144144       !-----------------------------------------------------------------------
    145        !   Computation of the soil temperatures using the Cgrd and Dgrd       
    146        !  coefficient computed at the previous time-step:                     
    147        !  -----------------------------------------------                     
    148        
    149        !    surface temperature                                               
    150        DO ig=1,ngrid 
     145       !   Computation of the soil temperatures using the Cgrd and Dgrd
     146       !  coefficient computed at the previous time-step:
     147       !  -----------------------------------------------
     148
     149       !    surface temperature
     150       DO ig=1,ngrid
    151151          ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/                   &
    152                &      (lambda*(1.-zd(ig,1))+1.)                                   
     152               &      (lambda*(1.-zd(ig,1))+1.)
    153153       ENDDO
    154        
    155        !   other temperatures                                                 
    156        DO jk=1,nsoil-1 
    157           DO ig=1,ngrid 
     154
     155       !   other temperatures
     156       DO jk=1,nsoil-1
     157          DO ig=1,ngrid
    158158             ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
    159159          ENDDO
    160160       ENDDO
    161        
     161
    162162    ENDIF
    163163
    164164    !-----------------------------------------------------------------------
    165     !   Computation of the Cgrd and Dgrd coefficient for the next step:     
    166     !   ---------------------------------------------------------------     
    167    
    168     DO jk=1,nsoil 
    169        zdz2(jk)=dz2(jk)/ptimestep 
    170     ENDDO
    171    
    172     DO ig=1,ngrid 
    173        z1(ig)=zdz2(nsoil)+dz1(nsoil-1) 
    174        zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig) 
    175        zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig) 
    176     ENDDO
    177    
    178     DO jk=nsoil-1,2,-1 
    179        DO ig=1,ngrid 
    180           z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk))) 
     165    !   Computation of the Cgrd and Dgrd coefficient for the next step:
     166    !   ---------------------------------------------------------------
     167
     168    DO jk=1,nsoil
     169       zdz2(jk)=dz2(jk)/ptimestep
     170    ENDDO
     171
     172    DO ig=1,ngrid
     173       z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
     174       zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
     175       zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
     176    ENDDO
     177
     178    DO jk=nsoil-1,2,-1
     179       DO ig=1,ngrid
     180          z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
    181181          zc(ig,jk-1)=                                                &
    182                &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)           
    183           zd(ig,jk-1)=dz1(jk-1)*z1(ig) 
     182               &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
     183          zd(ig,jk-1)=dz1(jk-1)*z1(ig)
    184184       ENDDO
    185185    ENDDO
    186    
    187     !-----------------------------------------------------------------------
    188     !   computation of the surface diffusive flux from ground and           
    189     !   calorific capacity of the ground:                                   
    190     !   ---------------------------------                                   
    191    
    192     DO ig=1,ngrid 
     186
     187    !-----------------------------------------------------------------------
     188    !   computation of the surface diffusive flux from ground and
     189    !   calorific capacity of the ground:
     190    !   ---------------------------------
     191
     192    DO ig=1,ngrid
    193193       pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*                              &
    194194            &   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
    … …  
    198198       pfluxgrd(ig)=pfluxgrd(ig)                                      &
    199199            &   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))   &
    200             &   /ptimestep                                                     
    201     ENDDO
    202    
     200            &   /ptimestep
     201    ENDDO
     202
    203203  END SUBROUTINE soil
    204  
     204
    205205END MODULE surface

  • dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90

    r4228 r4229  
    11MODULE turbulence
    2  
     2
    33#include "use_logging.h"
    44
    … …  
    88
    99  REAL, PARAMETER :: karman=0.4
    10   REAL :: lmixmin=100., emin_turb=1e-8 
     10  REAL :: lmixmin=100., emin_turb=1e-8
    1111
    1212  PUBLIC :: vdif, lmixmin, emin_turb
    13  
     13
    1414CONTAINS
    15  
     15
    1616  PURE SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
    1717    !=======================================================================
    … …  
    4747    !   Declarations:
    4848    !   -------------
    49    
     49
    5050    !   Arguments:
    5151    !   ----------
    52    
     52
    5353    INTEGER, INTENT(IN) ::  ngrid
    5454    REAL, INTENT(IN)  ::  pg, pz0(ngrid), pz(ngrid), &
    5555         pu(ngrid),pv(ngrid), pts(ngrid),ph(ngrid)
    5656    REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid)
    57    
     57
    5858    !   Local:
    5959    !   ------
    60    
     60
    6161    REAL, PARAMETER :: b=5., c=5., d=5., umin2=1e-12, &
    6262         c2b=2.*b, c3bc=3.*b*c, c3b=3.*b
    6363    INTEGER ig
    64     REAL zu2,z1,zri,zcd0,zz   
    65    
     64    REAL zu2,z1,zri,zcd0,zz
     65
    6666    !-----------------------------------------------------------------------
    6767    !   couche de surface:
    6868    !   ------------------
    69    
     69
    7070    !     DO ig=1,ngrid
    7171    !        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
    … …  
    7474    !     ENDDO
    7575    !     RETURN
    76    
    77    
     76
     77
    7878!!!! WARNING, verifier la formule originale de Louis!
    7979    DO ig=1,ngrid
    … …  
    9393       ENDIF
    9494    ENDDO
    95    
    96     !-----------------------------------------------------------------------
    97    
     95
     96    !-----------------------------------------------------------------------
     97
    9898  END SUBROUTINE vdif_cd
    99  
     99
    100100  PURE SUBROUTINE vdif_k(ngrid,nlay,   &
    101101       ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pkv,pkh)
    … …  
    106106         pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
    107107    REAL, INTENT(OUT) :: pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
    108    
     108
    109109    INTEGER ig,il
    110110    REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix
    111    
     111
    112112    DO ig=1,ngrid
    113113       pkv(ig,1)=0.
    … …  
    116116       pkh(ig,nlay+1)=0.
    117117    ENDDO
    118    
     118
    119119    DO il=2,nlay
    120120       DO ig=1,ngrid
    … …  
    138138       ENDDO
    139139    ENDDO
    140    
     140
    141141  END SUBROUTINE vdif_k
    142142
    143  
     143
    144144  SUBROUTINE multipl(n,x1,x2,y)
    145 !=======================================================================
    146 !   multiplication de deux vecteurs
    147 !=======================================================================
     145    !=======================================================================
     146    !   multiplication de deux vecteurs
     147    !=======================================================================
    148148    INTEGER n,i
    149149    REAL x1(n),x2(n),y(n)
    … …  
    162162       lwrite)
    163163    USE phys_const
    164    
    165 !=======================================================================
    166 !
    167 !   Diffusion verticale
    168 !   Shema implicite
    169 !   On commence par rajouter au variables x la tendance physique
    170 !   et on resoult en fait:
    171 !      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
    172 !
    173 !   arguments:
    174 !   ----------
    175 !
    176 !   entree:
    177 !   -------
    178 !
    179 !
    180 !=======================================================================
    181 
    182 !
    183 !   arguments:
    184 !   ----------
     164
     165    !=======================================================================
     166    !
     167    !   Diffusion verticale
     168    !   Shema implicite
     169    !   On commence par rajouter au variables x la tendance physique
     170    !   et on resoult en fait:
     171    !      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
     172    !
     173    !   arguments:
     174    !   ----------
     175    !
     176    !   entree:
     177    !   -------
     178    !
     179    !
     180    !=======================================================================
     181
     182    !
     183    !   arguments:
     184    !   ----------
    185185
    186186    INTEGER ngrid,nlay
    … …  
    198198    !   local:
    199199    !   ------
    200    
     200
    201201    INTEGER ilev,ig,ilay,nlev
    202202    INTEGER unit,ierr,it1,it2
    … …  
    213213    REAL zc(ngrid,nlay),zd(ngrid,nlay)
    214214    REAL zcst1
    215    
     215
    216216    !
    217217    !-----------------------------------------------------------------------
    218218    !   initialisations:
    219219    !   ----------------
    220    
     220
    221221    nlev=nlay+1
    222    
     222
    223223    !   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
    224224    !   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
    225225    !   ---------------------------------
    226    
     226
    227227    DO ilay=1,nlay
    228228       DO ig=1,ngrid
    … …  
    230230       ENDDO
    231231    ENDDO
    232    
     232
    233233    zcst1=4.*g*ptimestep/(r*r)
    234234    DO ilev=2,nlev-1
    … …  
    258258       LOG_DBG('vdif')
    259259    ENDIF
    260    
     260
    261261    !-----------------------------------------------------------------------
    262262    !   2. ajout des tendances physiques:
    263263    !   ------------------------------
    264    
     264
    265265    DO ilev=1,nlay
    266266       DO ig=1,ngrid
    … …  
    270270       ENDDO
    271271    ENDDO
    272    
     272
    273273    !-----------------------------------------------------------------------
    274274    !   3. calcul de  cd :
    … …  
    277277    CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
    278278    CALL vdif_k(ngrid,nlay,ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zkv,zkh)
    279    
     279
    280280    DO ig=1,ngrid
    281281       zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
    … …  
    283283       zcdh(ig)=zcdh(ig)*sqrt(zu2)
    284284    ENDDO
    285    
    286     IF(lwrite) THEN       
     285
     286    IF(lwrite) THEN
    287287       WRITELOG(*,*) 'Diagnostique diffusion verticale'
    288288       WRITELOG(*,*) 'LMIXMIN',lmixmin
    … …  
    295295       LOG_DBG('vdif')
    296296    ENDIF
    297    
     297
    298298    !-----------------------------------------------------------------------
    299299    !   integration verticale pour u:
    300300    !   -----------------------------
    301    
     301
    302302    CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
    303303    CALL multipl(ngrid,zcdv,zb0,zb)
    … …  
    307307       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
    308308    ENDDO
    309    
     309
    310310    DO ilay=nlay-1,1,-1
    311311       DO ig=1,ngrid
    … …  
    317317       ENDDO
    318318    ENDDO
    319    
     319
    320320    DO ig=1,ngrid
    321321       zu(ig,1)=zc(ig,1)
    … …  
    326326       ENDDO
    327327    ENDDO
    328    
     328
    329329    !-----------------------------------------------------------------------
    330330    !   integration verticale pour v:
    … …  
    336336       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
    337337    ENDDO
    338    
     338
    339339    DO ilay=nlay-1,1,-1
    340340       DO ig=1,ngrid
    … …  
    346346       ENDDO
    347347    ENDDO
    348    
     348
    349349    DO ig=1,ngrid
    350350       zv(ig,1)=zc(ig,1)
    … …  
    355355       ENDDO
    356356    ENDDO
    357    
     357
    358358    !-----------------------------------------------------------------------
    359359    !   integration verticale pour h:
    … …  
    367367       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
    368368    ENDDO
    369    
     369
    370370    DO ilay=nlay-1,1,-1
    371371       DO ig=1,ngrid
    … …  
    377377       ENDDO
    378378    ENDDO
    379    
     379
    380380    !-----------------------------------------------------------------------
    381381    !   rajout eventuel de planck dans le shema implicite:
    382382    !   --------------------------------------------------
    383    
     383
    384384    z4st=4.*5.67e-8*ptimestep
    385385    DO ig=1,ngrid
    386386       zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
    387387    ENDDO
    388    
     388
    389389    !-----------------------------------------------------------------------
    390390    !   calcul le l'evolution de la temperature du sol:
    391391    !   -----------------------------------------------
    392    
     392
    393393    DO ig=1,ngrid
    394394       z1(ig) = pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) &
    … …  
    399399       pdtsrf(ig) = (ztsrf2(ig)-ptsrf(ig))/ptimestep
    400400    ENDDO
    401    
     401
    402402    !-----------------------------------------------------------------------
    403403    !   integration verticale finale:
    404404    !   -----------------------------
    405    
     405
    406406    DO ilay=2,nlay
    407407       DO ig=1,ngrid
    … …  
    409409       ENDDO
    410410    ENDDO
    411    
     411
    412412    !-----------------------------------------------------------------------
    413413    !   calcul final des tendances de la diffusion verticale:
    414414    !   -----------------------------------------------------
    415    
     415
    416416    DO ilev = 1, nlay
    417417       DO ig=1,ngrid
    … …  
    424424       ENDDO
    425425    ENDDO
    426    
     426
    427427    IF(lwrite) THEN
    428        WRITELOG(*,*) 
     428       WRITELOG(*,*)
    429429       WRITELOG(*,*) 'Diagnostique de la diffusion verticale'
    430430       WRITELOG(*,*) 'h avant et apres diffusion verticale'
    … …  
    437437
    438438  END SUBROUTINE vdif
    439  
     439
    440440END MODULE turbulence

  • dynamico_lmdz/simple_physics/phyparam/physics/use_logging.h

    r4210 r4229  
    1    USE logging
     1  USE logging
    22#define WRITELOG(junk, fmt) logging_lineno = MIN(logging_bufsize, logging_lineno+1) ; WRITE(logging_buf(logging_lineno), fmt)
    33#define LOG_WARN(tag) CALL flush_log(log_level_warn, tag)

  • dynamico_lmdz/simple_physics/phyparam/physics/writefield_mod.F90

    r4227 r4229  
    2424  INTERFACE writefield
    2525     PROCEDURE writefield1, writefield2
    26   END INTERFACE writefield
     26  END INTERFACE
    2727
    2828  PUBLIC :: writefield
    … …  
    5353#endif
    5454  END SUBROUTINE writefield1
    55  
     55
    5656END MODULE writefield_mod
Note: See TracChangeset for help on using the changeset viewer.