Changeset 149 for codes


Ignore:
Timestamp:
04/03/13 12:05:12 (11 years ago)
Author:
sdubey
Message:
Added few new routines to read NC files and compute diagnostics to r145.
Few routines of dry physics including radiation module, surface process and convective adjustment in new routine phyparam.f90. dynetat to read start files for dynamics. check_conserve routine to compute conservation of quatities like mass, energy etc.etat0_heldsz.f90 for held-suarez test case initial conditions. new Key time_style=lmd or dcmip to use day_step, ndays like in LMDZ
Location:
codes/icosagcm/trunk/src
Files:
12 added
9 edited

Legend:

Unmodified
Added
Removed

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

    r148 r149  
    6666    CALL transfert_request(f_rhodz,req_i1) 
    6767         
    68     IF (is_mpi_root) PRINT *, 'Advection scheme' 
     68!    IF (is_mpi_root) PRINT *, 'Advection scheme' 
    6969 
    7070!    DO ind=1,ndomain 

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

    r146 r149  
    1111  TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:) 
    1212 
    13   PUBLIC init_caldyn, caldyn, write_output_fields 
     13  PUBLIC init_caldyn, caldyn, write_output_fields,un2ulonlat 
    1414 
    1515  INTEGER :: caldyn_hydrostat, caldyn_conserv 
    … …  
    722722     
    723723    CALL writefield("ps",f_ps) 
    724     CALL writefield("dps",f_dps) 
    725     CALL writefield("phis",f_phis) 
    726     CALL vorticity(f_u,f_buf_v) 
    727     CALL writefield("vort",f_buf_v) 
    728  
    729     CALL w_omega(f_ps, f_u, f_buf_i) 
    730     CALL writefield('omega', f_buf_i) 
    731     IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN 
    732       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level) 
    733       CALL writefield("omega"//TRIM(str_pression),f_buf_s) 
    734     ENDIF 
     724!    CALL writefield("dps",f_dps) 
     725!    CALL writefield("phis",f_phis) 
     726!    CALL vorticity(f_u,f_buf_v) 
     727!    CALL writefield("vort",f_buf_v) 
     728 
     729!    CALL w_omega(f_ps, f_u, f_buf_i) 
     730!    CALL writefield('omega', f_buf_i) 
     731!    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN 
     732!      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level) 
     733!      CALL writefield("omega"//TRIM(str_pression),f_buf_s) 
     734!    ENDIF 
    735735     
    736736    ! Temperature 
    … …  
    767767    ! geopotential 
    768768    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i) 
    769     CALL writefield("p",f_buf_p) 
    770     CALL writefield("phi",f_buf_i) 
    771     CALL writefield("theta",f_buf1_i) ! potential temperature 
    772     CALL writefield("pk",f_buf2_i) ! Exner pressure 
     769!    CALL writefield("p",f_buf_p) 
     770!    CALL writefield("phi",f_buf_i) 
     771     CALL writefield("theta",f_buf1_i) ! potential temperature 
     772!    CALL writefield("pk",f_buf2_i) ! Exner pressure 
    773773   
    774774   

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

    r148 r149  
    9696      IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2' 
    9797   CASE DEFAULT 
    98       IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip' 
     98      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), &  
     99        ' in dissip_gcm.f90/init_dissip' 
    99100      STOP 
    100101   END SELECT 

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

    r146 r149  
    271271            nf2=domain_glo(edge_glo(e)%assign_domain)%face 
    272272            d%edge_assign_sign(k,i,j)=1-2*MOD(12+tab_index(nf,nf2,0),2) 
    273             IF (MOD(6+k+tab_index(nf,nf2,0),6)/=edge_glo(e)%assign_pos .AND. MOD(6+k+tab_index(nf,nf2,0),6) /= MOD(edge_glo(e)%assign_pos+3,6)) THEN 
     273            IF (MOD(6+k+tab_index(nf,nf2,0),6)/=edge_glo(e)%assign_pos .AND. MOD(6+k+tab_index(nf,nf2,0),6) &  
     274                /= MOD(edge_glo(e)%assign_pos+3,6)) THEN 
    274275              d%edge_assign_sign(k,i,j)=-d%edge_assign_sign(k,i,j) 
    275276            ENDIF 

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

    r113 r149  
    11MODULE etat0_mod 
     2    CHARACTER(len=255),SAVE :: etat0_type 
    23 
    34CONTAINS 
    … …  
    1213    USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0   
    1314    USE etat0_dcmip5_mod, ONLY : etat0_dcmip5=>etat0   
     15    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0   
     16    USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0  
     17    USE dynetat0_hz_mod,  ONLY : dynetat0_hz=>etat0  
     18 
    1419    IMPLICIT NONE 
    1520    TYPE(t_field),POINTER :: f_ps(:) 
    … …  
    1924    TYPE(t_field),POINTER :: f_q(:) 
    2025     
    21     CHARACTER(len=255) :: etat0_type 
    2226    etat0_type='jablonowsky06' 
    2327    CALL getin("etat0",etat0_type) 
    … …  
    2832    CASE ('academic') 
    2933       CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     34    CASE ('heldsz') 
     35        print*,"heldsz test case" 
     36       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    3037    CASE ('dcmip1') 
    3138       CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    … …  
    3845     CASE ('dcmip5') 
    3946       CALL etat0_dcmip5(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     47     CASE ('readnf_start')  
     48          print*,"readnf_start used"     
     49       CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q)  
     50        CASE ('readnf_hz')  
     51          print*,"readnf_hz used" 
     52       CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q)  
    4053   CASE DEFAULT 
    4154       PRINT*, 'Bad selector for variable etat0 <',etat0_type, & 

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

    r131 r149  
    2525  CALL init_transfert 
    2626  CALL init_writefield 
     27 
    2728!  CALL allocate_field(sum_ne,field_T,type_real) 
    2829!  CALL allocate_field_glo(sum_ne_glo,field_T,type_real) 

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

    r99 r149  
    88  SUBROUTINE init_physics 
    99  USE icosa 
    10   USE physics_dcmip_mod, init_physics_dcmip=>init_physics 
     10  USE physics_dcmip_mod,init_physics_dcmip=>init_physics 
     11  USE physics_dry_mod 
    1112  IMPLICIT NONE 
    1213     
    … …  
    1819      CASE ('dcmip') 
    1920        CALL init_physics_dcmip 
     21 
     22      CASE ('lmd') 
     23        CALL init_physics_dry 
    2024       
    2125      CASE DEFAULT 
    22          PRINT*, 'Bad selector for variable physics <',physics_type, & 
     26         PRINT*, 'Bad selector for variable physics init <',physics_type, & 
    2327              '> options are <none>, <dcmip>,' 
    24          STOP 
     28 
    2529    END SELECT 
    2630     
    2731  END SUBROUTINE init_physics 
    2832 
    29   SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 
     33  SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 
    3034  USE icosa 
     35  USE physics_dry_mod 
    3136  USE physics_dcmip_mod, physics_dcmip=>physics 
     37  USE etat0_mod 
     38  USE etat0_heldsz_mod 
    3239  IMPLICIT NONE 
    3340    INTEGER, INTENT(IN)   :: it 
     41    REAL(rstd),INTENT(IN)::jD_cur,jH_cur 
    3442    TYPE(t_field),POINTER :: f_phis(:) 
    3543    TYPE(t_field),POINTER :: f_ps(:) 
    … …  
    3745    TYPE(t_field),POINTER :: f_ue(:) 
    3846    TYPE(t_field),POINTER :: f_q(:) 
     47    LOGICAL:: firstcall,lastcall 
    3948     
    4049    SELECT CASE(TRIM(physics_type)) 
    4150      CASE ('none') 
     51 
     52        SELECT CASE(TRIM(etat0_type)) 
     53        CASE('heldsz')  
     54!       CALL transfert_request(f_ps,req_i1) 
     55!       CALL transfert_request(f_theta_rhodz,req_i1) 
     56!       CALL transfert_request(f_ue,req_e1_vect) 
     57!       CALL held_saurez(f_ps,f_theta_rhodz,f_ue)  
     58        CASE DEFAULT 
     59        PRINT*,"NO PHYSICAL PACAKAGE USED"  
     60        END SELECT 
    4261     
    4362      CASE ('dcmip') 
    4463        CALL physics_dcmip(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 
     64 
     65      CASE ('dry') 
     66        CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 
    4567       
    4668      CASE DEFAULT 
    4769         PRINT*, 'Bad selector for variable physics <',physics_type, & 
    4870              '> options are <none>, <dcmip>,' 
    49          STOP 
     71  STOP 
    5072    END SELECT 
    5173     

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

    r132 r149  
    1111  INTEGER,SAVE    :: itau_out, itau_adv, itau_dissip, itau_physics, itaumax 
    1212   
     13  INTEGER,SAVE :: day_step,ndays 
     14  REAL(rstd),SAVE :: jD_ref,jH_ref 
     15  INTEGER,SAVE :: day_ini,day_end,annee_ref,day_ref 
     16  REAL(rstd),SAVE::start_time 
     17  CHARACTER(LEN=255) :: time_style 
     18  INTEGER,SAVE:: an, mois, jour 
     19  REAL(rstd),SAVE:: heure 
     20  CHARACTER (LEN=10):: calend 
     21 
    1322  PUBLIC create_time_counter_header, update_time_counter, close_time_counter, init_time,  & 
    14          dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax 
     23         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax, &  
     24day_step,ndays,jD_ref,jH_ref,day_ini,day_end,annee_ref,day_ref,an, mois, jour,heure, & 
     25            calend,time_style 
     26 
    1527 
    1628 
    … …  
    2436  IMPLICIT NONE 
    2537  REAL(rstd) :: run_length 
    26    
     38 
     39 
     40     time_style='dcmip' 
     41   CALL getin('time_style',time_style) 
     42 
     43   IF (TRIM(time_style)=='dcmip')  Then 
    2744    dt=90. 
    2845    CALL getin('dt',dt) 
    … …  
    3047    itaumax=100 
    3148    CALL getin('itaumax',itaumax) 
    32  
    33     itau_adv=1 
    34     CALL getin('itau_adv',itau_adv) 
    35  
    36     itau_dissip=1 
    37     CALL getin('itau_dissip',itau_dissip) 
    38  
    39     itau_physics=1 
    40     CALL getin('itau_physics',itau_physics) 
    4149 
    4250    run_length=dt*itaumax 
    … …  
    5462    itau_out=FLOOR(.5+write_period/dt) 
    5563    IF (is_mpi_root) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 
     64    ENDIF  
     65 
     66    itau_adv=1 
     67    CALL getin('itau_adv',itau_adv) 
     68 
     69    itau_dissip=1 
     70    CALL getin('itau_dissip',itau_dissip) 
     71 
     72    itau_physics=1 
     73    CALL getin('itau_physics',itau_physics) 
     74 
    5675     
    5776    CALL create_time_counter_header 

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

    r148 r149  
    1919  USE physics_mod 
    2020  USE mpipara 
     21  USE IOIPSL 
     22  USE maxicosa 
     23  USE check_conserve_mod 
    2124  USE trace 
    2225  USE transfert_mod 
    … …  
    4346  REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:) 
    4447  REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 
    45  
    4648!  REAL(rstd) :: dt, run_length 
    4749  INTEGER :: ind 
    … …  
    4951  CHARACTER(len=255) :: scheme_name 
    5052  LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 
     53  CHARACTER(len=7) :: first  
     54  REAL(rstd),SAVE :: jD_cur, jH_cur 
     55  REAL(rstd),SAVE :: start_time 
    5156  LOGICAL, PARAMETER :: check=.FALSE. 
    5257!  INTEGER :: itaumax 
    … …  
    8893  CALL allocate_field(f_wfluxt,field_t,type_real,llm+1) 
    8994 
     95!---------------------------------------------------- 
     96  IF (TRIM(time_style)=='lmd')  Then 
     97 
     98   day_step=180 
     99   CALL getin('day_step',day_step) 
     100 
     101   ndays=1 
     102   CALL getin('ndays',ndays) 
     103 
     104   dt = daysec/REAL(day_step) 
     105   itaumax = ndays*day_step 
     106 
     107   calend = 'earth_360d' 
     108   CALL getin('calend', calend) 
     109 
     110   day_ini = 0 
     111   CALL getin('day_ini',day_ini) 
     112 
     113   day_end = 0 
     114   CALL getin('day_end',day_end) 
     115 
     116   annee_ref = 1998 
     117   CALL getin('annee_ref',annee_ref) 
     118 
     119   start_time = 0 
     120   CALL getin('start_time',start_time)  
     121 
     122   write_period=0 
     123   CALL getin('write_period',write_period) 
     124       
     125   write_period=write_period/scale_factor 
     126   itau_out=FLOOR(write_period/dt) 
     127    
     128   PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out  
     129 
     130  mois = 1 ; heure = 0. 
     131  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref) 
     132  jH_ref = jD_ref - int(jD_ref)  
     133  jD_ref = int(jD_ref)  
     134   
     135  CALL ioconf_startdate(INT(jD_ref),jH_ref)  
     136  write(*,*)'annee_ref, mois, day_ref, heure, jD_ref' 
     137  write(*,*)annee_ref, mois, day_ref, heure, jD_ref 
     138  write(*,*)"ndays,day_step,itaumax,dt======>" 
     139  write(*,*)ndays,day_step,itaumax,dt 
     140  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure) 
     141  write(*,*)'jD_ref+jH_ref,an, mois, jour, heure' 
     142  write(*,*)jD_ref+jH_ref,an, mois, jour, heure 
     143  day_end = day_ini + ndays  
     144 END IF  
     145!---------------------------------------------------- 
     146 
    90147  scheme_name='runge_kutta' 
    91148  CALL getin('scheme',scheme_name) 
    … …  
    126183!  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm) 
    127184!  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm) 
    128  
    129185!  CALL allocate_field(f_theta,field_t,type_real,llm) 
    130186!  CALL allocate_field(f_dtheta,field_t,type_real,llm) 
    … …  
    135191  CALL init_advect_tracer 
    136192  CALL init_physics 
    137  
    138     
    139   CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     193!========================================= INITIALIZATION  
     194! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
     195  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
    140196  CALL writefield("phis",f_phis,once=.TRUE.) 
    141197  CALL transfert_request(f_q,req_i1)  
    … …  
    158214     CALL compute_rhodz(.FALSE., ps, rhodz)    
    159215  END DO 
    160    
     216 
    161217  CALL transfert_request(f_phis,req_i0)  
     218  CALL transfert_request(f_phis,req_i1)  
    162219  CALL transfert_request(f_phis,req_i1)  
    163220 
    … …  
    170227    ENDIF 
    171228     
    172     IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It 
    173229    IF (mod(it,itau_out)==0 ) THEN 
    174       CALL writefield("q",f_q) 
     230!      IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It 
    175231      CALL update_time_counter(dt*it) 
     232      CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)   
    176233    ENDIF 
    177      
    178     CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
     234 
     235      CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
    179236 
    180237    DO stage=1,nb_stage 
    … …  
    189246       CASE (mlf) 
    190247          CALL  leapfrog_matsuno_scheme(stage) 
    191            
    192248          !      CASE ('leapfrog') 
    193           !        CALL leapfrog_scheme 
     249          !      CALL leapfrog_scheme 
    194250          ! 
    195251          !      CASE ('adam_bashforth') 
    196           !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
    197           !        CALL adam_bashforth_scheme 
     252          !      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) 
     253          !      CALL adam_bashforth_scheme 
    198254       CASE DEFAULT 
    199255          STOP 
    … …  
    208264    IF(MOD(it+1,itau_adv)==0) THEN 
    209265!       CALL transfert_request(f_wfluxt,req_i1) ! FIXME 
    210 !       CALL transfert_request(f_hfluxt,req_e1) ! FIXME 
     266!      CALL transfert_request(f_hfluxt,req_e1) ! FIXME 
    211267 
    212268       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step 
    … …  
    227283         END DO 
    228284       ENDIF 
    229  
    230285    END IF 
    231  
     286!---------------------------------------------------- 
     287  jD_cur = jD_ref + day_ini - day_ref + it/day_step 
     288  jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step) 
     289  jD_cur = jD_cur + int(jH_cur) 
     290  jH_cur = jH_cur - int(jH_cur) 
     291!       print*,"Just b4 phys" 
     292    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q) 
     293!---------------------------------------------------- 
    232294!    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q) 
    233295    ENDDO 
    … …  
    493555          DO i=ii_begin-dd,ii_end+dd 
    494556             ij=(j-1)*iim+i 
    495              m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g  
     557             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
    496558             IF(comp) THEN 
    497559                rhodz(ij,l) = m 
    … …  
    508570          STOP 
    509571       ELSE 
    510 !          PRINT *, 'No discrepancy between ps and rhodz detected' 
     572          PRINT *, 'No discrepancy between ps and rhodz detected' 
    511573       END IF 
    512574    END IF 
Note: See TracChangeset for help on using the changeset viewer.