Changeset 346


Ignore:
Timestamp:
08/01/15 01:52:51 (9 years ago)
Author:
dubos
Message:

Cleanup DCMIP4

Location:
codes/icosagcm/trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/etat0.f90

    r345 r346  
    1919    USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0 
    2020    USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0 
     21    USE etat0_dcmip4_mod, ONLY : getin_etat0_dcmip4=>getin_etat0 
    2122    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0 
    2223    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0 
    … …  
    2425    ! Ad hoc interfaces 
    2526    USE etat0_academic_mod, ONLY : etat0_academic=>etat0   
    26     USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0   
    2727    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0   
    2828    USE etat0_venus_mod,  ONLY : etat0_venus=>etat0   
    … …  
    6565       CALL getin_etat0_dcmip2 
    6666    CASE ('dcmip3') 
     67    CASE ('dcmip4') 
     68        CALL getin_etat0_dcmip4 
    6769    CASE ('dcmip5') 
    6870        CALL getin_etat0_dcmip5 
    … …  
    8688       CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q) 
    8789       PRINT *, "Venus (Lebonnois et al., 2012) test case" 
    88     CASE ('dcmip4') 
    89         IF(nqtot<2) THEN 
    90            IF (is_mpi_root)  THEN 
    91               PRINT *, "nqtot must be at least 2 for test case DCMIP4" 
    92            END IF 
    93            STOP 
    94         END IF 
    95         CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    9690   CASE DEFAULT 
    9791      IF(collocated) THEN 
    … …  
    170164    USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0 
    171165    USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0 
     166    USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0 
    172167    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0 
    173168    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0 
    … …  
    213208       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    214209       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
     210    CASE('dcmip4') 
     211       CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
     212       CALL compute_dcmip4(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    215213    CASE('dcmip5') 
    216214       CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 

  • codes/icosagcm/trunk/src/etat0_dcmip4.f90

    r286 r346  
    11MODULE etat0_dcmip4_mod 
    22  USE icosa 
     3  IMPLICIT NONE 
    34  PRIVATE 
     5   
    46  REAL(rstd),PARAMETER :: eta0=0.252 
    57  REAL(rstd),PARAMETER :: etat=0.2 
    … …  
    1113  REAL(rstd),PARAMETER :: Gamma=0.005 
    1214  REAL(rstd),PARAMETER :: up0=1 
    13   REAL(rstd) :: latw=2*Pi/9 
    14 !$OMP THREADPRIVATE(latw) 
    15   REAL(rstd) :: pw=34000 
    16 !$OMP THREADPRIVATE(pw) 
    17   REAL(rstd) :: q0=0.021 
    18 !$OMP THREADPRIVATE(q0) 
    19   REAL(rstd) :: lonc 
    20 !$OMP THREADPRIVATE(lonc) 
    21   REAL(rstd) :: latc 
    22 !$OMP THREADPRIVATE(latc) 
     15  REAL(rstd),PARAMETER :: lonc=Pi/9, latc=2*Pi/9, latw=2*Pi/9 
     16  REAL(rstd),PARAMETER :: pw=34000 
     17  REAL(rstd),PARAMETER :: q0=0.021 
     18   
     19  INTEGER,SAVE :: testcase 
     20  !$OMP THREADPRIVATE(testcase)   
     21   
     22  PUBLIC getin_etat0, compute_etat0 
    2323 
    24   INTEGER,SAVE :: testcase 
    25 !$OMP THREADPRIVATE(testcase)   
     24CONTAINS 
    2625 
    27   PUBLIC  etat0 
    28 CONTAINS 
    29    
    30     
    31      
    32   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    33   USE icosa 
    34   IMPLICIT NONE 
    35     TYPE(t_field),POINTER :: f_ps(:) 
    36     TYPE(t_field),POINTER :: f_phis(:) 
    37     TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    38     TYPE(t_field),POINTER :: f_u(:) 
    39     TYPE(t_field),POINTER :: f_q(:) 
    40    
    41     REAL(rstd),POINTER :: ps(:) 
    42     REAL(rstd),POINTER :: phis(:) 
    43     REAL(rstd),POINTER :: theta_rhodz(:,:) 
    44     REAL(rstd),POINTER :: u(:,:) 
    45     REAL(rstd),POINTER :: q(:,:,:) 
    46     INTEGER :: ind 
    47  
     26  SUBROUTINE getin_etat0 
     27    USE mpipara, ONLY : is_mpi_root 
     28    IF(nqtot<2) THEN 
     29       IF (is_mpi_root)  THEN 
     30          PRINT *, "nqtot must be at least 2 for test case DCMIP4" 
     31       END IF 
     32       STOP 
     33    END IF 
    4834    testcase=1 
    4935    CALL getin("dcmip4_testcase",testcase) 
     36  END SUBROUTINE getin_etat0 
     37 
     38  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q) 
     39    USE icosa 
     40    USE disvert_mod 
     41    INTEGER, INTENT(IN) :: ngrid 
     42    REAL(rstd),INTENT(IN) :: lon(ngrid) 
     43    REAL(rstd),INTENT(IN) :: lat(ngrid) 
     44    REAL(rstd),INTENT(OUT) :: phis(ngrid) 
     45    REAL(rstd),INTENT(OUT) :: ps(ngrid) 
     46    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm) 
     47    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) 
     48    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) 
     49    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) 
    5050     
    51     DO ind=1,ndomain 
    52       IF (.NOT. assigned_domain(ind)) CYCLE 
    53       CALL swap_dimensions(ind) 
    54       CALL swap_geometry(ind) 
    55       ps=f_ps(ind) 
    56       phis=f_phis(ind) 
    57       theta_rhodz=f_theta_rhodz(ind) 
    58       u=f_u(ind) 
    59       q=f_q(ind) 
    60       CALL compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q) 
    61     ENDDO 
    62  
    63   END SUBROUTINE etat0 
    64    
    65   SUBROUTINE compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q) 
    66   USE icosa 
    67   USE disvert_mod 
    68   USE pression_mod 
    69   USE exner_mod 
    70   USE geopotential_mod 
    71   USE theta2theta_rhodz_mod 
    72   IMPLICIT NONE   
    73   REAL(rstd),INTENT(OUT) :: ps(iim*jjm) 
    74   REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 
    75   REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 
    76   REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 
    77   REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) 
    78    
    79   INTEGER :: i,j,l,ij 
    80   REAL(rstd) :: theta(iim*jjm,llm) 
    81   REAL(rstd) :: Y(iim*jjm,llm) 
    82   REAL(rstd) :: vort 
    83   REAL(rstd) :: eta(llm) 
    84   REAL(rstd) :: etav(llm) 
    85   REAL(rstd) :: etas, etavs 
    86   REAL(rstd) :: lon,lat 
    87   REAL(rstd) :: ulon(3) 
    88   REAL(rstd) :: ep(3), norm_ep 
    89   REAL(rstd) :: Tave, T 
    90   REAL(rstd) :: phis_ave 
    91   REAL(rstd) :: V0(3) 
    92   REAL(rstd) :: r2 
    93   REAL(rstd) :: utot 
    94   REAL(rstd) :: lonx,latx 
    95   REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r 
    96    
    97     lonc=Pi/9 
    98     latc=2*Pi/9  
    99    
    100     DO l=1,llm 
    101       eta(l)= 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) ) 
    102       etav(l)=(eta(l)-eta0)*Pi/2 
    103     ENDDO 
     51    INTEGER :: l,ij 
     52    REAL(rstd) :: etal, etavl, etas, etavs, sinlat, coslat, & 
     53         Y, Tave, T, phis_ave, vort, r2, utot, & 
     54         dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r 
     55     
    10456    etas=ap(1)/preff+bp(1) 
    10557    etavs=(etas-eta0)*Pi/2 
    106  
    107     DO j=jj_begin,jj_end 
    108       DO i=ii_begin,ii_end 
    109         ij=(j-1)*iim+i 
    110         ps(ij)=ps0 
    111       ENDDO 
     58    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) 
     59     
     60    DO ij=1,ngrid 
     61       sinlat=SIN(lat(ij)) 
     62       coslat=COS(lat(ij)) 
     63       phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sinlat**6 * (coslat**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  & 
     64            +(8./5*coslat**3 * (sinlat**2 + 2./3) - Pi/4)*radius*Omega ) 
     65       ps(ij)=ps0 
    11266    ENDDO 
    11367     
     68    DO l=1,llm 
     69       etal = 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) ) 
     70       etavl=(etal-eta0)*Pi/2 
     71        
     72       Tave=T0*etal**(Rd*Gamma/g) 
     73       dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* etal**(Rd*Gamma/g-kappa-1) 
     74       IF (etat>etal) THEN 
     75          Tave=Tave+DeltaT*(etat-etal)**5 
     76          dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-etal)**4 * etal**(-kappa)  & 
     77               + kappa * (etat-etal)**5 * etal**(-kappa-1)) 
     78       END IF 
     79        
     80       DO ij=1,ngrid 
     81          sinlat=SIN(lat(ij)) 
     82          coslat=COS(lat(ij)) 
     83           
     84          K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc) 
     85          r=radius*acos(K) 
     86          utot=u0*cos(etavl)**1.5*sin(2*lat(ij))**2 + up0*exp(-(r/(0.1*radius))**2) 
     87          ulon(ij,l) = utot 
     88          ulat(ij,l) = 0. 
     89          Y = ((-2*sinlat**6*(coslat**2+1./3)+10./63)*2*u0*cos(etavl)**1.5     & 
     90               + (8./5*coslat**3*(sinlat**2+2./3)-Pi/4)*radius*Omega) 
     91          T = Tave + 0.75*(etal*Pi*u0/Rd)*sin(etavl)*cos(etavl)**0.5 * Y 
     92          temp(ij,l)=T 
     93           
     94          IF (testcase==1) THEN 
     95             q(ij,l,1)=T*etal**(-kappa) 
     96             IF(nqtot>2) q(ij,l,3)=1. 
     97             dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*etal**(-kappa)*sin(etavl)*cos(etavl)**0.5 * Y &  
     98                  + 3/8. * Pi**2*u0/Rd * etal**(1-kappa) * cos(etavl)**1.5 * Y                    &  
     99                  - 3./16. * Pi**2 * u0 /Rd * etal**(1-kappa) * sin(etavl)**2 * cos(etavl)**(-0.5) *Y & 
     100                  - 9./8.  * Pi**2 * u0 /Rd * etal**(1-kappa) * sin(etavl)**2 * cos(etavl)   & 
     101                  * (-2*sinlat**6*(coslat**2+1./3.)+10./63.)  
     102             dthetaodlat=3./4.*Pi*u0/Rd*etal**(1-kappa)*sin(etavl)*cos(etavl)**0.5                                   & 
     103                  *( 2*u0*cos(etavl)**1.5 * ( -12 * coslat*sinlat**5*(coslat**2+1./3.)+4*coslat*sinlat**7) & 
     104                  + radius*omega*(-24./5. * sinlat * coslat**2 * (sinlat**2 + 2./3.) + 16./5. * coslat**4 * sinlat)) 
     105              
     106             duodeta=-u0 * sin(2*lat(ij))**2 * 3./4.*Pi * cos(etavl)**0.5 * sin(etavl) 
     107              
     108             vort = -4*u0/radius*cos(etavl)**1.5 * sinlat * coslat * (2.-5.*sinlat**2)                  &  
     109                  + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat(ij))-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*coslat & 
     110                  -cos(latc)*sinlat*cos(lon(ij)-lonc))/(sqrt(1-K**2))) 
     111             q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sinlat*omega+vort)*dthetaodeta)) 
     112             IF(nqtot>3) q(ij,l,4)=cos(lon(ij))*coslat 
     113          ELSE IF (testcase==2) THEN 
     114             q(ij,l,1)=q0*exp(-(lat(ij)/latw)**4)*exp(-((etal-1)*preff/pw)**2) 
     115          END IF 
     116       END DO 
     117    END DO 
    114118     
    115     CALL lonlat2xyz(lonc,latc,V0) 
    116  
    117     u(:,:)=1e10       
    118     DO l=1,llm 
    119       DO j=jj_begin-1,jj_end+1 
    120         DO i=ii_begin-1,ii_end+1 
    121           ij=(j-1)*iim+i 
    122            
    123           lon=lon_e(ij+u_right) ; lat=lat_e(ij+u_right) 
    124           K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) 
    125           r=radius*acos(K) 
    126           utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) 
    127           u(ij+u_right,l) = utot * sum(elon_e(ij+u_right,:) * ep_e(ij+u_right,:)) 
    128  
    129           lon=lon_e(ij+u_lup) ; lat=lat_e(ij+u_lup) 
    130           K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) 
    131           r=radius*acos(K) 
    132           utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) 
    133           u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:)) 
    134  
    135           lon=lon_e(ij+u_ldown) ; lat=lat_e(ij+u_ldown) 
    136           K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) 
    137           r=radius*acos(K) 
    138           utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2) 
    139           u(ij+u_ldown,l) = utot * sum(elon_e(ij+u_ldown,:) * ep_e(ij+u_ldown,:)) 
    140  
    141        ENDDO 
    142      ENDDO 
    143     ENDDO  
    144       
    145       
    146      DO l=1,llm 
    147        Tave=T0*eta(l)**(Rd*Gamma/g) 
    148        IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 
    149        DO j=jj_begin-1,jj_end+1 
    150          DO i=ii_begin-1,ii_end+1 
    151            ij=(j-1)*iim+i 
    152            lat=lat_i(ij) 
    153             Y(ij,l)=((-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5     & 
    154                                 + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega) 
    155             T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) 
    156              
    157             theta(ij,l)=T*eta(l)**(-kappa) 
    158  
    159           ENDDO 
    160        ENDDO 
    161      ENDDO 
    162       
    163       
    164      phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) 
    165      DO j=jj_begin,jj_end 
    166        DO i=ii_begin,ii_end 
    167          ij=(j-1)*iim+i 
    168          lat=lat_i(ij) 
    169          phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  & 
    170                                            +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega ) 
    171 !         phis(ij)=phis_ave+u0*cos(etavs)**1.5 
    172  
    173        ENDDO 
    174      ENDDO 
    175  
    176     CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) 
    177  
    178  
    179     IF (testcase==1) THEN 
    180      
    181       q(:,:,1)=theta(:,:) 
    182       IF(nqtot>2) q(:,:,3)=1. 
    183  
    184       DO l=1,llm 
    185          dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1) 
    186          IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa)  & 
    187                                                           + kappa * (etat-eta(l))**5 * eta(l)**(-kappa-1)) 
    188          DO j=jj_begin,jj_end 
    189            DO i=ii_begin,ii_end 
    190              ij=(j-1)*iim+i 
    191              lon=lon_i(ij) ; lat=lat_i(ij) 
    192              dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) &  
    193                                         + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l)                    &  
    194                                    - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)& 
    195                                    - 9./8.  * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))   & 
    196                                                                                   * (-2*sin(lat)**6*(cos(lat)**2+1./3.)+10./63.)  
    197              dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5                                   & 
    198                          *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) & 
    199                         + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat))) 
    200                            
    201              duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l)) 
    202             
    203              K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) 
    204              r=radius*acos(K) 
    205              vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2)                  &  
    206                     + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) & 
    207                                                             -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2))) 
    208              q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta)) 
    209              IF(nqtot>3) q(ij,l,4)=cos(lon)*cos(lat) 
    210            ENDDO 
    211          ENDDO 
    212        ENDDO        
    213  
    214     ELSE IF (testcase==2) THEN 
    215      
    216       DO l=1,llm 
    217          DO j=jj_begin,jj_end 
    218            DO i=ii_begin,ii_end 
    219              ij=(j-1)*iim+i 
    220              q(ij,l,1)=q0*exp(-(lat_i(ij)/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2) 
    221            ENDDO 
    222          ENDDO 
    223        ENDDO 
    224      
    225     ENDIF        
    226  
    227          
    228   END SUBROUTINE compute_etat0_dcmip4 
     119  END SUBROUTINE compute_etat0 
    229120   
    230121END MODULE etat0_dcmip4_mod 
Note: See TracChangeset for help on using the changeset viewer.