Changeset 344


Ignore:
Timestamp:
07/31/15 17:29:55 (9 years ago)
Author:
dubos
Message:

Cleanup DCMIP 1&2

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

Legend:

Unmodified
Added
Removed

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

    r327 r344  
    11MODULE etat0_mod 
    22  USE icosa 
     3  IMPLICIT NONE          
    34  PRIVATE 
    45 
    … …  
    1617    USE disvert_mod 
    1718    ! New interface 
     19    USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0 
     20    USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0 
    1821    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0 
    1922    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0 
    … …  
    2124    ! Old interface 
    2225    USE etat0_academic_mod, ONLY : etat0_academic=>etat0   
    23     USE etat0_dcmip1_mod, ONLY : etat0_dcmip1=>etat0 
    24     USE etat0_dcmip2_mod, ONLY : etat0_dcmip2=>etat0 
    2526    USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0   
    2627    USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0   
    … …  
    3839     
    3940    REAL(rstd),POINTER :: ps(:), mass(:,:) 
    40     LOGICAL :: init_mass 
     41    LOGICAL :: init_mass, collocated 
    4142    INTEGER :: ind,i,j,ij,l 
    4243 
    … …  
    5253    CALL getin("etat0",etat0_type) 
    5354     
     55    !------------------- New interface --------------------- 
     56    collocated=.TRUE. 
    5457    SELECT CASE (TRIM(etat0_type)) 
    55        !------------------- New interface --------------------- 
    5658    CASE ('isothermal') 
    5759       CALL getin_etat0_isothermal 
    58        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 
    5960    CASE ('temperature_profile') 
    6061       CALL getin_etat0_temperature 
    61        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 
    6262    CASE ('jablonowsky06') 
    63        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 
    64      CASE ('dcmip5') 
     63    CASE ('dcmip1') 
     64        CALL getin_etat0_dcmip1 
     65    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') 
     66       CALL getin_etat0_dcmip2 
     67    CASE ('dcmip5') 
    6568        CALL getin_etat0_dcmip5 
    66         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 
    6769    CASE ('williamson91.6') 
    6870       init_mass=.FALSE. 
    6971       CALL getin_etat0_williamson 
    70        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 
    71        !------------------- Old interface -------------------- 
     72    CASE DEFAULT 
     73       collocated=.FALSE. 
     74    END SELECT 
     75 
     76    !------------------- Old interface -------------------- 
     77    SELECT CASE (TRIM(etat0_type)) 
    7278    CASE ('start_file') 
    7379       CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    … …  
    8086       CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q) 
    8187       PRINT *, "Venus (Lebonnois et al., 2012) test case" 
    82     CASE ('dcmip1') 
    83        CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    84     CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') 
    85        CALL etat0_dcmip2(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    8688    CASE ('dcmip3') 
    8789       CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    … …  
    9597        CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    9698   CASE DEFAULT 
    97        PRINT*, 'Bad selector for variable etat0 <',etat0_type, & 
    98             '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> ' 
    99        STOP 
     99      IF(collocated) THEN 
     100         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 
     101      ELSE 
     102         PRINT*, 'Bad selector for variable etat0 <',etat0_type, & 
     103              '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> ' 
     104         STOP 
     105      END IF 
    100106    END SELECT 
    101107 
    … …  
    163169    USE wind_mod 
    164170    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0 
     171    USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0 
     172    USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0 
    165173    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0 
    166174    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0 
    … …  
    197205       CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i) 
    198206       CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e) 
     207    CASE('dcmip1') 
     208       CALL compute_dcmip1(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
     209       CALL compute_dcmip1(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
     210    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') 
     211       CALL compute_dcmip2(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i) 
     212       CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)       
    199213    CASE('dcmip5') 
    200214       CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 

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

    r286 r344  
    11MODULE etat0_dcmip1_mod 
    22  USE icosa 
     3  IMPLICIT NONE          
    34  PRIVATE 
     5 
     6  INTEGER, PARAMETER :: const=1, cos_bell=2, slotted_cyl=3, & 
     7       dbl_cos_bell_q1=4, dbl_cos_bell_q2=5, complement=6, hadley=7, & 
     8       dcmip11=-1 
     9  INTEGER, SAVE :: shape 
     10!$OMP THREADPRIVATE(shape) 
    411 
    512  REAL(rstd), SAVE  :: h0=1. 
    … …  
    3643!$OMP THREADPRIVATE(zc) 
    3744 
    38   PUBLIC etat0 
     45  PUBLIC getin_etat0, compute_etat0 
    3946   
    4047CONTAINS 
    41    
    42   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
    43   USE icosa 
    44   IMPLICIT NONE 
    45     TYPE(t_field),POINTER :: f_ps(:) 
    46     TYPE(t_field),POINTER :: f_phis(:) 
    47     TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    48     TYPE(t_field),POINTER :: f_u(:) 
    49     TYPE(t_field),POINTER :: f_q(:) 
    50    
    51     REAL(rstd),POINTER :: ps(:) 
    52     REAL(rstd),POINTER :: phis(:) 
    53     REAL(rstd),POINTER :: theta_rhodz(:,:) 
    54     REAL(rstd),POINTER :: u(:,:) 
    55     REAL(rstd),POINTER :: q(:,:,:) 
     48 
     49  SUBROUTINE getin_etat0 
    5650    CHARACTER(len=255) :: dcmip1_adv_shape 
    57     INTEGER :: ind 
    58      
    5951    R0=radius*0.5 
    6052    rt=radius*0.5 
    6153    dcmip1_adv_shape='cos_bell' 
    6254    CALL getin('dcmip1_shape',dcmip1_adv_shape) 
    63      
    64     DO ind=1,ndomain 
    65       IF (.NOT. assigned_domain(ind)) CYCLE 
    66       CALL swap_dimensions(ind) 
    67       CALL swap_geometry(ind) 
    68       ps=f_ps(ind) 
    69       phis=f_phis(ind) 
    70       theta_rhodz=f_theta_rhodz(ind) 
    71       u=f_u(ind) 
    72       q=f_q(ind) 
    73  
    74  
    75       SELECT CASE(TRIM(dcmip1_adv_shape)) 
    76       CASE('const') 
    77          CALL compute_etat0_ncar(1,ps, phis, theta_rhodz, u, q(:,:,1)) 
    78       CASE('cos_bell') 
    79          CALL compute_etat0_ncar(2,ps, phis, theta_rhodz, u, q(:,:,1)) 
    80       CASE('slotted_cyl') 
    81          CALL compute_etat0_ncar(3,ps, phis, theta_rhodz, u, q(:,:,1)) 
    82       CASE('dbl_cos_bell_q1') 
    83          CALL compute_etat0_ncar(4,ps, phis, theta_rhodz, u, q(:,:,1)) 
    84       CASE('dbl_cos_bell_q2') 
    85          CALL compute_etat0_ncar(5,ps, phis, theta_rhodz, u, q(:,:,1)) 
    86       CASE('complement') 
    87          CALL compute_etat0_ncar(6,ps, phis, theta_rhodz, u, q(:,:,1)) 
    88       CASE('hadley')  ! hadley like meridional circulation  
    89          CALL compute_etat0_ncar(7,ps, phis, theta_rhodz, u, q(:,:,1)) 
    90       CASE('dcmip11') 
    91          IF(nqtot==5) THEN 
    92             CALL compute_etat0_ncar(4,ps, phis, theta_rhodz, u, q(:,:,1)) 
    93             CALL compute_etat0_ncar(5,ps, phis, theta_rhodz, u, q(:,:,2)) 
    94             CALL compute_etat0_ncar(3,ps, phis, theta_rhodz, u, q(:,:,3)) 
    95             CALL compute_etat0_ncar(6,ps, phis, theta_rhodz, u, q(:,:,4)) 
    96             CALL compute_etat0_ncar(1,ps, phis, theta_rhodz, u, q(:,:,5)) 
    97          ELSE 
    98             PRINT *,'Error : etat0_dcmip=dcmip11 and nqtot = ',nqtot,' .' 
    99             PRINT *,'nqtot must be equal to 5 when etat0_dcmip=dcmip11' 
    100             STOP 
    101          END IF 
    102       CASE DEFAULT 
    103          PRINT *, 'Bad selector for variable dcmip1_adv_shape : <', TRIM(dcmip1_adv_shape),   & 
    104               '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', & 
    105               '<dbl_cos_bell_q2>, <complement>, <hadley>' 
    106          STOP 
    107       END SELECT 
    108  
    109     ENDDO 
    110  
    111   END SUBROUTINE etat0 
    112    
    113   SUBROUTINE compute_etat0_ncar(icase, ps, phis, theta_rhodz, u, q) 
     55    SELECT CASE(TRIM(dcmip1_adv_shape)) 
     56    CASE('const') 
     57       shape=const 
     58    CASE('cos_bell') 
     59       shape=cos_bell 
     60    CASE('slotted_cyl') 
     61       shape=slotted_cyl 
     62    CASE('dbl_cos_bell_q1') 
     63       shape=dbl_cos_bell_q1 
     64    CASE('dbl_cos_bell_q2') 
     65       shape=dbl_cos_bell_q2 
     66    CASE('complement') 
     67       shape=complement 
     68    CASE('hadley')  ! hadley like meridional circulation  
     69       shape=hadley 
     70    CASE('dcmip11') 
     71       IF(nqtot<5) THEN 
     72          PRINT *,'Error : etat0_dcmip=dcmip11 and nqtot = ',nqtot,' .' 
     73          PRINT *,'nqtot must be equal to 5 when etat0_dcmip=dcmip11' 
     74          STOP 
     75       END IF 
     76       shape=dcmip11 
     77    CASE DEFAULT 
     78       PRINT *, 'Bad selector for variable dcmip1_adv_shape : <', TRIM(dcmip1_adv_shape),   & 
     79            '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', & 
     80            '<dbl_cos_bell_q2>, <complement>, <hadley>' 
     81       STOP 
     82    END SELECT 
     83  END SUBROUTINE getin_etat0 
     84 
     85  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q) 
     86    INTEGER, INTENT(IN) :: ngrid 
     87    REAL(rstd),INTENT(IN) :: lon(ngrid),lat(ngrid) 
     88    REAL(rstd),INTENT(OUT) :: ps(ngrid),phis(ngrid) 
     89    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm),ulon(ngrid,llm),ulat(ngrid,llm) 
     90    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) 
     91    ps = ncar_p0 
     92    phis = 0. ; temp = 0. ;  
     93    ulon = 0. ; ulat=0. 
     94    SELECT CASE(shape) 
     95    CASE(dcmip11) 
     96       CALL compute_etat0_ncar(4,ngrid,lon,lat,q(:,:,1)) 
     97       CALL compute_etat0_ncar(5,ngrid,lon,lat,q(:,:,2)) 
     98       CALL compute_etat0_ncar(3,ngrid,lon,lat,q(:,:,3)) 
     99       CALL compute_etat0_ncar(6,ngrid,lon,lat,q(:,:,4)) 
     100       CALL compute_etat0_ncar(1,ngrid,lon,lat,q(:,:,5)) 
     101    CASE DEFAULT 
     102       CALL compute_etat0_ncar(shape,ngrid,lon,lat,q(:,:,1)) 
     103    END SELECT 
     104  END SUBROUTINE compute_etat0 
     105 
     106  SUBROUTINE compute_etat0_ncar(icase,ngrid,lon,lat, q) 
    114107  USE icosa 
    115108  USE disvert_mod 
    116   USE pression_mod 
    117   USE exner_mod 
    118   USE geopotential_mod 
    119   USE theta2theta_rhodz_mod 
    120109  IMPLICIT NONE   
    121   INTEGER, INTENT(in) :: icase 
    122   REAL(rstd),INTENT(OUT) :: ps(iim*jjm) 
    123   REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 
    124   REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 
    125   REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 
    126   REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm)   
    127  
    128   REAL(rstd) :: qxt1(iim*jjm,llm) 
    129   REAL(rstd) :: lon, lat 
    130   REAL(rstd) ::dd1,dd2,dd1t1,dd1t2,dd2t1 
    131   REAL(rstd) :: pr, zr(llm+1), zrl(llm) 
    132   REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx 
    133   REAL(rstd) :: X2(3),X1(3) 
    134   INTEGER :: i,j,n,l 
    135  
    136   u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0 
     110  INTEGER, INTENT(IN) :: icase, ngrid 
     111  REAL(rstd),INTENT(IN) :: lon(ngrid),lat(ngrid)   
     112  REAL(rstd),INTENT(OUT) :: q(ngrid,llm)   
     113  REAL(rstd) :: zr(llm+1), zrl(llm), qxt1(ngrid,llm) 
     114  REAL(rstd) :: pr 
     115  !  REAL(rstd) :: lon, lat 
     116  INTEGER :: n,l 
    137117   
    138118  DO l=1, llm+1 
    … …  
    149129     q=1 
    150130  CASE(2) 
    151      !--------------------------------------------- SINGLE COSINE BELL 
    152131     CALL cosine_bell_1(q)       
    153132  CASE(3) 
    154133     CALL slotted_cylinders(q)   
    155134  CASE(4) 
    156      PRINT *, 'Double cosine bell' 
    157135     CALL cosine_bell_2(q)       
    158136  CASE(5) 
    … …  
    176154  END SELECT 
    177155       
    178 CONTAINS 
    179  
     156  CONTAINS 
    180157!====================================================================== 
    181158 
    182159    SUBROUTINE cosine_bell_1(hx) 
    183160    IMPLICIT NONE  
    184     REAL(rstd) :: hx(iim*jjm,llm) 
     161    REAL(rstd) :: hx(ngrid,llm) 
     162    REAL(rstd) :: rr1,rr2    
     163    INTEGER :: n,l 
     164    DO l=1,llm  
     165       DO n=1,ngrid 
     166          CALL dist_lonlat(lon0,lat0,lon(n),lat(n),rr1)  ! GC distance from center  
     167          rr1 = radius*rr1 
     168          IF ( rr1 .LT. R0 ) then  
     169             hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0))  
     170          ELSE 
     171             hx(n,l)=0.0 
     172          END IF 
     173       END DO 
     174    END DO 
     175  END SUBROUTINE cosine_bell_1 
     176 
     177!============================================================================== 
     178   
     179  SUBROUTINE cosine_bell_2(hx)  
     180    REAL(rstd) :: hx(ngrid,llm) 
     181    REAL(rstd) :: rr1,rr2,dd1,dd2,dd1t1,dd1t2,dd2t1 
     182    INTEGER :: n,l 
     183    DO l=1,llm 
     184       DO n=1,ngrid 
     185          CALL dist_lonlat(lonc1,latc1,lon(n),lat(n),rr1)  ! GC distance from center  
     186          rr1 = radius*rr1 
     187          CALL dist_lonlat(lonc2,latc2,lon(n),lat(n),rr2)  ! GC distance from center  
     188          rr2 = radius*rr2 
     189          dd1t1 = rr1/rt 
     190          dd1t1 = dd1t1*dd1t1 
     191          dd1t2 = (zrl(l) - zc)/zt 
     192          dd1t2 = dd1t2*dd1t2  
     193          dd1 = dd1t1 + dd1t2  
     194          dd1 = Min(1.0,dd1)   
     195          dd2t1 = rr2/rt 
     196          dd2t1 = dd2t1*dd2t1 
     197          dd2 = dd2t1 + dd1t2  
     198          dd2 = Min(1.0,dd2)  
     199          hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2)) 
     200       END DO 
     201    END DO   
     202  END SUBROUTINE cosine_bell_2 
     203   
     204  !============================================================================= 
     205 
     206  SUBROUTINE slotted_cylinders(hx)  
     207    REAL(rstd) :: hx(ngrid,llm) 
     208    REAL(rstd) :: rr1,rr2,dd1,dd2,dd1t1,dd1t2,dd2t1 
     209    INTEGER :: n,l 
     210    DO l=1,llm 
     211       DO n=1,ngrid 
     212          CALL dist_lonlat(lonc1,latc1,lon(n),lat(n),rr1)  ! GC distance from center  
     213          rr1 = radius*rr1 
     214          CALL dist_lonlat(lonc2,latc2,lon(n),lat(n),rr2)  ! GC distance from center  
     215          rr2 = radius*rr2 
     216          dd1t1 = rr1/rt 
     217          dd1t1 = dd1t1*dd1t1 
     218          dd1t2 = (zrl(l) - zc)/zt 
     219          dd1t2 = dd1t2*dd1t2  
     220          dd1 = dd1t1 + dd1t2   
     221          dd2t1 = rr2/rt 
     222          dd2t1 = dd2t1*dd2t1 
     223          dd2 = dd2t1 + dd1t2 
     224          IF ( dd1 .LT. 0.5 ) Then  
     225             hx(n,l) = 1.0 
     226          ELSEIF ( dd2 .LT. 0.5 ) Then  
     227             hx(n,l) = 1.0 
     228          ELSE 
     229             hx(n,l) = 0.1  
     230          END IF 
     231          IF ( zrl(l) .GT. zc ) Then  
     232             IF ( ABS(latc1 - lat_i(n)) .LT. 0.125 ) Then  
     233                hx(n,l)= 0.1 
     234             ENDIF 
     235          ENDIF 
     236       END DO 
     237    END DO 
     238  END SUBROUTINE slotted_cylinders 
    185239     
    186       DO l=1,llm  
    187         DO j=jj_begin-1,jj_end+1 
    188           DO i=ii_begin-1,ii_end+1 
    189             n=(j-1)*iim+i 
    190             CALL dist_lonlat(lon0,lat0,lon_i(n),lat_i(n),rr1)  ! GC distance from center  
    191             rr1 = radius*rr1 
    192             IF ( rr1 .LT. R0 ) then  
    193               hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0))  
    194             ELSE 
    195               hx(n,l)=0.0 
    196             END IF 
    197           END DO 
    198         END DO 
    199       END DO 
    200     
    201     END SUBROUTINE cosine_bell_1 
    202  
    203240!============================================================================== 
    204241 
    205     SUBROUTINE   cosine_bell_2(hx)  
    206     IMPLICIT NONE        
    207     REAL(rstd) :: hx(iim*jjm,llm) 
    208     
    209       DO l=1,llm 
    210         DO j=jj_begin-1,jj_end+1 
    211           DO i=ii_begin-1,ii_end+1 
    212             n=(j-1)*iim+i 
    213             CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center  
    214             rr1 = radius*rr1 
    215             CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center  
    216             rr2 = radius*rr2 
    217             dd1t1 = rr1/rt 
    218             dd1t1 = dd1t1*dd1t1 
    219             dd1t2 = (zrl(l) - zc)/zt 
    220             dd1t2 = dd1t2*dd1t2  
    221             dd1 = dd1t1 + dd1t2  
    222             dd1 = Min(1.0,dd1)   
    223             dd2t1 = rr2/rt 
    224             dd2t1 = dd2t1*dd2t1 
    225             dd2 = dd2t1 + dd1t2  
    226             dd2 = Min(1.0,dd2)  
    227             hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2)) 
    228           END DO 
    229         END DO 
    230       END DO 
    231  
    232     END SUBROUTINE cosine_bell_2 
    233  
    234 !============================================================================= 
    235  
    236     SUBROUTINE slotted_cylinders(hx)  
    237     IMPLICIT NONE 
    238     REAL(rstd) :: hx(iim*jjm,llm) 
    239      
    240       DO l=1,llm 
    241         DO j=jj_begin-1,jj_end+1 
    242           DO i=ii_begin-1,ii_end+1 
    243             n=(j-1)*iim+i 
    244             CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center  
    245             rr1 = radius*rr1 
    246             CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center  
    247             rr2 = radius*rr2 
    248             dd1t1 = rr1/rt 
    249             dd1t1 = dd1t1*dd1t1 
    250             dd1t2 = (zrl(l) - zc)/zt 
    251             dd1t2 = dd1t2*dd1t2  
    252             dd1 = dd1t1 + dd1t2   
    253             dd2t1 = rr2/rt 
    254             dd2t1 = dd2t1*dd2t1 
    255             dd2 = dd2t1 + dd1t2 
    256              
    257             IF ( dd1 .LT. 0.5 ) Then  
    258               hx(n,l) = 1.0 
    259             ELSEIF ( dd2 .LT. 0.5 ) Then  
    260               hx(n,l) = 1.0 
    261             ELSE 
    262               hx(n,l) = 0.1  
    263             END IF   
    264      
    265             IF ( zrl(l) .GT. zc ) Then  
    266               IF ( ABS(latc1 - lat) .LT. 0.125 ) Then  
    267                 hx(n,l)= 0.1 
    268               ENDIF  
    269             ENDIF   
    270    
    271           ENDDO 
    272        END DO 
    273       END DO 
    274  
    275     END SUBROUTINE slotted_cylinders  
    276  
    277 !============================================================================== 
    278  
    279242    SUBROUTINE hadleyq(hx)  
    280       IMPLICIT NONE  
    281       REAL(rstd)::hx(iim*jjm,llm)  
     243      REAL(rstd)::hx(ngrid,llm)  
    282244      REAL(rstd),PARAMETER:: zz1=2000.,zz2=5000.,zz0=0.5*(zz1+zz2) 
     245      INTEGER :: n,l 
    283246       
    284247      DO l=1,llm 

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

    r295 r344  
    33! test cases DCMIP 2012, category 2 : Orographic gravity waves 
    44 
    5   USE genmod 
    6  
     5  USE icosa 
    76  PRIVATE 
    87 
    9   PUBLIC  etat0 
     8  INTEGER, SAVE :: testcase 
     9  INTEGER, PARAMETER :: mountain=0, schaer_noshear=1, schaer_shear=2 
     10 
     11  PUBLIC getin_etat0, compute_etat0 
    1012 
    1113CONTAINS 
    1214 
    13  
    14   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    15   USE icosa 
    16   USE theta2theta_rhodz_mod 
    17   IMPLICIT NONE 
    18     TYPE(t_field),POINTER :: f_ps(:) 
    19     TYPE(t_field),POINTER :: f_phis(:) 
    20     TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    21     TYPE(t_field),POINTER :: f_u(:) 
    22     TYPE(t_field),POINTER :: f_q(:) 
    23          
    24     TYPE(t_field),POINTER,SAVE :: f_temp(:) 
    25     REAL(rstd),POINTER :: ps(:) 
    26     REAL(rstd),POINTER :: phis(:) 
    27     REAL(rstd),POINTER :: u(:,:) 
    28     REAL(rstd),POINTER :: Temp(:,:) 
    29     REAL(rstd),POINTER :: q(:,:,:) 
    30  
    31     INTEGER :: ind, icase 
     15  SUBROUTINE getin_etat0 
    3216    CHARACTER(len=255) :: etat0_type 
    3317    etat0_type='jablonowsky06' 
    3418    CALL getin("etat0",etat0_type) 
    35  
    3619    SELECT CASE (TRIM(etat0_type)) 
    3720    CASE('dcmip2_mountain') 
    38        icase=0 
     21       testcase = mountain 
    3922    CASE('dcmip2_schaer_noshear') 
    40        icase=1 
     23       testcase = schaer_noshear 
    4124    CASE('dcmip2_schaer_shear') 
    42        icase=2 
     25       testcase = schaer_shear 
    4326    CASE DEFAULT 
    4427       PRINT *, 'This should not happen : etat0_type =', TRIM(etat0_type), ' in etat0_dcmip2.f90/etat0' 
    … …  
    4629    END SELECT 
    4730    PRINT *, 'Orographic gravity-wave test case :', TRIM(etat0_type) 
     31  END SUBROUTINE getin_etat0 
    4832 
    49     CALL allocate_field(f_temp,field_t,type_real,llm,name='temp')  
     33  SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, Temp, ulon, ulat) 
     34    USE disvert_mod 
     35    IMPLICIT NONE 
     36    INTEGER, INTENT(IN)    :: ngrid 
     37    REAL(rstd),INTENT(IN)  :: lon(ngrid) 
     38    REAL(rstd),INTENT(IN)  :: lat(ngrid) 
     39    REAL(rstd),INTENT(OUT) :: ps(ngrid) 
     40    REAL(rstd),INTENT(OUT) :: phis(ngrid) 
     41    REAL(rstd),INTENT(OUT) :: Temp(ngrid,llm) 
     42    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) 
     43    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) 
     44    INTEGER :: l,ij 
     45    REAL(rstd) :: hyam, hybm 
    5046 
    51     DO ind=1,ndomain 
    52       IF (.NOT. assigned_domain(ind)) CYCLE 
    53       CALL swap_dimensions(ind) 
    54       CALL swap_geometry(ind) 
    55  
    56       ps=f_ps(ind) 
    57       phis=f_phis(ind) 
    58       u=f_u(ind) 
    59       q=f_q(ind) 
    60       temp=f_temp(ind) 
    61       CALL compute_etat0_DCMIP2(icase,ps,phis,u,Temp,q) 
    62     ENDDO 
    63      
    64     CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) 
    65              
    66     CALL deallocate_field(f_temp) 
    67     
    68   END SUBROUTINE etat0 
    69    
    70   SUBROUTINE compute_etat0_DCMIP2(icase, ps, phis, u, temp,q) 
    71     USE icosa 
    72     USE disvert_mod, ONLY : ap,bp 
    73     USE pression_mod 
    74     USE theta2theta_rhodz_mod 
    75     USE wind_mod 
    76     IMPLICIT NONE 
    77     INTEGER, INTENT(IN) :: icase 
    78     REAL(rstd), INTENT(OUT) :: ps(iim*jjm) 
    79     REAL(rstd), INTENT(OUT) :: phis(iim*jjm) 
    80     REAL(rstd), INTENT(OUT) :: u(3*iim*jjm,llm) 
    81     REAL(rstd), INTENT(OUT) :: temp(iim*jjm,llm) 
    82     REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm) 
    83     REAL(rstd) :: ulon(3*iim*jjm,llm), ulat(3*iim*jjm,llm) 
    84      
    85     INTEGER :: i,j,l,ij 
    86     REAL(rstd) :: hyam, hybm, pp, dummy1, dummy2, dummy3   
    87      
    8847    ! Hexagons : ps,phis,temp 
    8948    DO l=1,llm 
    … …  
    9150       hyam = .5*(ap(l)+ap(l+1))/preff 
    9251       hybm = .5*(bp(l)+bp(l+1)) 
    93        DO j=jj_begin,jj_end 
    94           DO i=ii_begin,ii_end 
    95              ij=(j-1)*iim+i 
    96              CALL comp_all(lon_i(ij),lat_i(ij), ps(ij),phis(ij),temp(ij,l), dummy1,dummy2) 
    97           END DO 
     52       DO ij=1,ngrid 
     53          CALL comp_all(hyam,hybm,lon(ij),lat(ij), & 
     54               ps(ij),phis(ij), Temp(ij,l),ulon(ij,l), ulat(ij,l) ) 
    9855       END DO 
    9956    END DO 
    100      
    101     ! Edges : ulon,ulat 
    102     DO l=1,llm 
    103        ! The surface pressure is not set yet so we provide the hybrid coefficients 
    104        hyam = .5*(ap(l)+ap(l+1))/preff 
    105        hybm = .5*(bp(l)+bp(l+1)) 
    106        DO j=jj_begin,jj_end 
    107           DO i=ii_begin,ii_end 
    108              ij=(j-1)*iim+i 
    109              CALL comp_all(lon_e(ij+u_right), lat_e(ij+u_right), & 
    110                   dummy1,dummy2,dummy3, ulon(ij+u_right,l),ulat(ij+u_right,l)) 
    111              CALL comp_all(lon_e(ij+u_lup), lat_e(ij+u_lup), & 
    112                   dummy1,dummy2,dummy3, ulon(ij+u_lup,l),ulat(ij+u_lup,l)) 
    113              CALL comp_all(lon_e(ij+u_ldown), lat_e(ij+u_ldown), & 
    114                   dummy1,dummy2,dummy3, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l)) 
    115           END DO 
    116        END DO 
    117     END DO 
    118     CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u) 
     57  END SUBROUTINE compute_etat0 
    11958 
    120     q=1. 
    121  
    122     CONTAINS 
    123        
    124       SUBROUTINE comp_all(lon,lat, psj,phisj,tempj, ulonj,ulatj) 
    125         USE dcmip_initial_conditions_test_1_2_3 
    126         REAL(rstd), INTENT(IN) :: lon, lat 
    127         REAL(rstd), INTENT(OUT) :: psj,phisj,tempj,ulonj,ulatj 
    128         REAL :: dummy 
    129         dummy=0. 
    130         SELECT CASE (icase) 
    131         CASE(0) 
    132            CALL test2_steady_state_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm, & 
    133                 ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy) 
    134         CASE(1) 
    135            CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,0,& 
    136                 ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy) 
    137         CASE(2) 
    138            CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,1, & 
    139                 ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy) 
    140         END SELECT 
    141       END SUBROUTINE comp_all 
    142  
    143   END SUBROUTINE compute_etat0_DCMIP2 
    144  
     59  SUBROUTINE comp_all(hyam,hybm,lon,lat, psj,phisj,tempj, ulonj,ulatj) 
     60    USE dcmip_initial_conditions_test_1_2_3 
     61    REAL(rstd), INTENT(IN) :: hyam, hybm, lon, lat 
     62    REAL(rstd), INTENT(OUT) :: psj,phisj,tempj,ulonj,ulatj 
     63    REAL :: dummy 
     64    dummy=0. 
     65    SELECT CASE (testcase) 
     66    CASE(mountain) 
     67       CALL test2_steady_state_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm, & 
     68            ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy) 
     69    CASE(schaer_noshear) 
     70       CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,0,& 
     71            ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy) 
     72    CASE(schaer_shear) 
     73       CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,1, & 
     74            ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy) 
     75    END SELECT 
     76  END SUBROUTINE comp_all 
    14577 
    14678END MODULE etat0_dcmip2_mod 
Note: See TracChangeset for help on using the changeset viewer.