Index: time.f90 =================================================================== --- time.f90 (revision 142) +++ time.f90 (working copy) @@ -10,11 +10,23 @@ REAL(rstd),SAVE :: write_period INTEGER,SAVE :: itau_out, itau_adv, itau_dissip, itau_physics, itaumax + INTEGER,SAVE :: day_step,ndays + REAL(rstd),SAVE :: jD_ref,jH_ref + INTEGER,SAVE :: day_ini,day_end,annee_ref,day_ref + REAL(rstd),SAVE::start_time + CHARACTER(LEN=255) :: time_style + INTEGER,SAVE:: an, mois, jour + REAL(rstd),SAVE:: heure + CHARACTER (LEN=10):: calend + PUBLIC create_time_counter_header, update_time_counter, close_time_counter, init_time, & - dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax + dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax, & +day_step,ndays,jD_ref,jH_ref,day_ini,day_end,annee_ref,day_ref,an, mois, jour,heure, & + calend,time_style + CONTAINS SUBROUTINE init_time @@ -23,22 +35,18 @@ USE mpipara IMPLICIT NONE REAL(rstd) :: run_length - + + + time_style='dcmip' + CALL getin('time_style',time_style) + + IF (TRIM(time_style)=='dcmip') Then dt=90. CALL getin('dt',dt) itaumax=100 CALL getin('itaumax',itaumax) - itau_adv=1 - CALL getin('itau_adv',itau_adv) - - itau_dissip=1 - CALL getin('itau_dissip',itau_dissip) - - itau_physics=1 - CALL getin('itau_physics',itau_physics) - run_length=dt*itaumax CALL getin('run_length',run_length) itaumax=run_length/dt @@ -53,6 +61,17 @@ write_period=write_period/scale_factor itau_out=FLOOR(.5+write_period/dt) IF (is_mpi_root) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out + ENDIF + + itau_adv=1 + CALL getin('itau_adv',itau_adv) + + itau_dissip=1 + CALL getin('itau_dissip',itau_dissip) + + itau_physics=1 + CALL getin('itau_physics',itau_physics) + CALL create_time_counter_header Index: caldyn_gcm.f90 =================================================================== --- caldyn_gcm.f90 (revision 142) +++ caldyn_gcm.f90 (working copy) @@ -10,7 +10,7 @@ TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:) TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:) - PUBLIC init_caldyn, caldyn, write_output_fields + PUBLIC init_caldyn, caldyn, write_output_fields,un2ulonlat INTEGER :: caldyn_hydrostat, caldyn_conserv @@ -708,17 +708,17 @@ str_pression=ADJUSTL(str_pression) CALL writefield("ps",f_ps) - CALL writefield("dps",f_dps) - CALL writefield("phis",f_phis) - CALL vorticity(f_u,f_buf_v) - CALL writefield("vort",f_buf_v) +! CALL writefield("dps",f_dps) +! CALL writefield("phis",f_phis) +! CALL vorticity(f_u,f_buf_v) +! CALL writefield("vort",f_buf_v) - CALL w_omega(f_ps, f_u, f_buf_i) - CALL writefield('omega', f_buf_i) - IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN - CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level) - CALL writefield("omega"//TRIM(str_pression),f_buf_s) - ENDIF +! CALL w_omega(f_ps, f_u, f_buf_i) +! CALL writefield('omega', f_buf_i) +! IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN +! CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level) +! CALL writefield("omega"//TRIM(str_pression),f_buf_s) +! ENDIF ! Temperature CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; @@ -753,10 +753,10 @@ ! geopotential CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i) - CALL writefield("p",f_buf_p) - CALL writefield("phi",f_buf_i) - CALL writefield("theta",f_buf1_i) ! potential temperature - CALL writefield("pk",f_buf2_i) ! Exner pressure +! CALL writefield("p",f_buf_p) +! CALL writefield("phi",f_buf_i) + CALL writefield("theta",f_buf1_i) ! potential temperature +! CALL writefield("pk",f_buf2_i) ! Exner pressure END SUBROUTINE write_output_fields Index: advect_tracer.f90 =================================================================== --- advect_tracer.f90 (revision 142) +++ advect_tracer.f90 (working copy) @@ -57,7 +57,7 @@ CALL transfert_request(f_q,req_i1) CALL transfert_request(f_rhodz,req_i1) - IF (is_mpi_root) PRINT *, 'Advection scheme' +! IF (is_mpi_root) PRINT *, 'Advection scheme' ! DO ind=1,ndomain ! CALL swap_dimensions(ind) Index: dissip_gcm.f90 =================================================================== --- dissip_gcm.f90 (revision 142) +++ dissip_gcm.f90 (working copy) @@ -93,7 +93,8 @@ rayleigh_friction_type=2 IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2' CASE DEFAULT - IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip' + IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & + ' in dissip_gcm.f90/init_dissip' STOP END SELECT Index: physics.f90 =================================================================== --- physics.f90 (revision 142) +++ physics.f90 (working copy) @@ -7,7 +7,8 @@ SUBROUTINE init_physics USE icosa - USE physics_dcmip_mod, init_physics_dcmip=>init_physics + USE physics_dcmip_mod,init_physics_dcmip=>init_physics + USE physics_dry_mod IMPLICIT NONE CALL getin("physics",physics_type) @@ -17,36 +18,57 @@ CASE ('dcmip') CALL init_physics_dcmip + + CASE ('lmd') + CALL init_physics_dry CASE DEFAULT - PRINT*, 'Bad selector for variable physics <',physics_type, & + PRINT*, 'Bad selector for variable physics init <',physics_type, & '> options are , ,' - STOP + END SELECT END SUBROUTINE init_physics - SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_q) + SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q) USE icosa + USE physics_dry_mod USE physics_dcmip_mod, physics_dcmip=>physics + USE etat0_mod + USE etat0_heldsz_mod IMPLICIT NONE INTEGER, INTENT(IN) :: it + REAL(rstd),INTENT(IN)::jD_cur,jH_cur TYPE(t_field),POINTER :: f_phis(:) TYPE(t_field),POINTER :: f_ps(:) TYPE(t_field),POINTER :: f_theta_rhodz(:) TYPE(t_field),POINTER :: f_ue(:) TYPE(t_field),POINTER :: f_q(:) + LOGICAL:: firstcall,lastcall SELECT CASE(TRIM(physics_type)) CASE ('none') + + SELECT CASE(TRIM(etat0_type)) + CASE('heldsz') + CALL transfert_request(f_ps,req_i1) + CALL transfert_request(f_theta_rhodz,req_i1) + CALL transfert_request(f_ue,req_e1) + CALL held_saurez(f_ps,f_theta_rhodz,f_ue) + CASE DEFAULT + PRINT*,"NO PHYSICAL PACAKAGE USED" + END SELECT CASE ('dcmip') CALL physics_dcmip(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) + + CASE ('dry') + CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q) CASE DEFAULT PRINT*, 'Bad selector for variable physics <',physics_type, & '> options are , ,' - STOP + STOP END SELECT END SUBROUTINE physics Index: etat0.f90 =================================================================== --- etat0.f90 (revision 142) +++ etat0.f90 (working copy) @@ -1,4 +1,5 @@ MODULE etat0_mod + CHARACTER(len=255),SAVE :: etat0_type CONTAINS @@ -11,6 +12,10 @@ USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0 USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0 USE etat0_dcmip5_mod, ONLY : etat0_dcmip5=>etat0 + USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0 + USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0 + USE dynetat0_hz_mod, ONLY : dynetat0_hz=>etat0 + IMPLICIT NONE TYPE(t_field),POINTER :: f_ps(:) TYPE(t_field),POINTER :: f_phis(:) @@ -18,7 +23,6 @@ TYPE(t_field),POINTER :: f_u(:) TYPE(t_field),POINTER :: f_q(:) - CHARACTER(len=255) :: etat0_type etat0_type='jablonowsky06' CALL getin("etat0",etat0_type) @@ -27,6 +31,9 @@ CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CASE ('academic') CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q) + CASE ('heldsz') + print*,"heldsz test case" + CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CASE ('dcmip1') CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') @@ -37,6 +44,12 @@ CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CASE ('dcmip5') CALL etat0_dcmip5(f_ps,f_phis,f_theta_rhodz,f_u, f_q) + CASE ('readnf_start') + print*,"readnf_start used" + CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q) + CASE ('readnf_hz') + print*,"readnf_hz used" + CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q) CASE DEFAULT PRINT*, 'Bad selector for variable etat0 <',etat0_type, & '> options are , , ' Index: timeloop_gcm.f90 =================================================================== --- timeloop_gcm.f90 (revision 142) +++ timeloop_gcm.f90 (working copy) @@ -17,7 +17,11 @@ USE advect_tracer_mod USE physics_mod USE mpipara - +!------------------------------- + USE IOIPSL + USE maxicosa + USE check_conserve_mod +!------------------------------ IMPLICIT NONE TYPE(t_field),POINTER :: f_phis(:) ! TYPE(t_field),POINTER :: f_theta(:) @@ -46,6 +50,10 @@ INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme CHARACTER(len=255) :: scheme_name LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time + CHARACTER(len=7) :: first + REAL(rstd),SAVE :: jD_cur, jH_cur + REAL(rstd),SAVE :: start_time + ! INTEGER :: itaumax ! REAL(rstd) ::write_period ! INTEGER :: itau_out @@ -84,6 +92,58 @@ CALL allocate_field(f_hfluxt,field_u,type_real,llm) ! mass "fluxes" accumulated in time CALL allocate_field(f_wfluxt,field_t,type_real,llm+1) +!---------------------------------------------------- + IF (TRIM(time_style)=='lmd') Then + + day_step=180 + CALL getin('day_step',day_step) + + ndays=1 + CALL getin('ndays',ndays) + + dt = daysec/REAL(day_step) + itaumax = ndays*day_step + + calend = 'earth_360d' + CALL getin('calend', calend) + + day_ini = 0 + CALL getin('day_ini',day_ini) + + day_end = 0 + CALL getin('day_end',day_end) + + annee_ref = 1998 + CALL getin('annee_ref',annee_ref) + + start_time = 0 + CALL getin('start_time',start_time) + + write_period=0 + CALL getin('write_period',write_period) + + write_period=write_period/scale_factor + itau_out=FLOOR(write_period/dt) + + PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out + + mois = 1 ; heure = 0. + call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref) + jH_ref = jD_ref - int(jD_ref) + jD_ref = int(jD_ref) + + CALL ioconf_startdate(INT(jD_ref),jH_ref) + write(*,*)'annee_ref, mois, day_ref, heure, jD_ref' + write(*,*)annee_ref, mois, day_ref, heure, jD_ref + write(*,*)"ndays,day_step,itaumax,dt======>" + write(*,*)ndays,day_step,itaumax,dt + call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure) + write(*,*)'jD_ref+jH_ref,an, mois, jour, heure' + write(*,*)jD_ref+jH_ref,an, mois, jour, heure + day_end = day_ini + ndays + END IF +!---------------------------------------------------- + scheme_name='runge_kutta' CALL getin('scheme',scheme_name) @@ -122,7 +182,6 @@ ! CALL allocate_field(f_dum2,field_u,type_real,llm) ! CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm) ! CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm) - ! CALL allocate_field(f_theta,field_t,type_real,llm) ! CALL allocate_field(f_dtheta,field_t,type_real,llm) @@ -132,7 +191,11 @@ CALL init_advect_tracer CALL init_physics - CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) +!============================ SET 0 to all field + CALL initial0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) +! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) + CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) +!================================================ CALL writefield("phis",f_phis,once=.TRUE.) CALL transfert_request(f_q,req_i1) @@ -153,17 +216,18 @@ rhodz=f_rhodz(ind); ps=f_ps(ind) CALL compute_rhodz(.FALSE., ps, rhodz) END DO + DO it=0,itaumax - IF (is_mpi_root) PRINT *,"It No :",It," t :",dt*It IF (mod(it,itau_out)==0 ) THEN - CALL writefield("q",f_q) +! IF (is_mpi_root) PRINT *,"It No :",It," t :",dt*It CALL update_time_counter(dt*it) + CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) ENDIF - - CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) + CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) + DO stage=1,nb_stage CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & f_phis,f_ps,f_theta_rhodz,f_u, f_q, & @@ -177,11 +241,11 @@ CALL leapfrog_matsuno_scheme(stage) ! CASE ('leapfrog') - ! CALL leapfrog_scheme + ! CALL leapfrog_scheme ! ! CASE ('adam_bashforth') - ! CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) - ! CALL adam_bashforth_scheme + ! CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz) + ! CALL adam_bashforth_scheme CASE DEFAULT STOP END SELECT @@ -192,7 +256,7 @@ IF(MOD(it+1,itau_adv)==0) THEN CALL transfert_request(f_wfluxt,req_i1) ! FIXME -! CALL transfert_request(f_hfluxt,req_e1) ! FIXME +! CALL transfert_request(f_hfluxt,req_e1) ! FIXME CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step fluxt_zero=.TRUE. @@ -207,11 +271,17 @@ CALL swap_geometry(ind) rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) - CALL compute_rhodz(.FALSE., ps, rhodz) +!!!!! CALL compute_rhodz(.FALSE., ps, rhodz) END DO - END IF - +!---------------------------------------------------- + jD_cur = jD_ref + day_ini - day_ref + it/day_step + jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step) + jD_cur = jD_cur + int(jH_cur) + jH_cur = jH_cur - int(jH_cur) +! print*,"Just b4 phys" + CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q) +!---------------------------------------------------- ! CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q) ENDDO @@ -399,7 +469,7 @@ DO j=jj_begin-dd,jj_end+dd DO i=ii_begin-dd,ii_end+dd ij=(j-1)*iim+i - m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g + m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g IF(comp) THEN rhodz(ij,l) = m ELSE @@ -414,7 +484,7 @@ PRINT *, 'Discrepancy between ps and rhodz detected', err STOP ELSE -! PRINT *, 'No discrepancy between ps and rhodz detected' + PRINT *, 'No discrepancy between ps and rhodz detected' END IF END IF Index: icosa_gcm.f90 =================================================================== --- icosa_gcm.f90 (revision 142) +++ icosa_gcm.f90 (working copy) @@ -24,6 +24,7 @@ CALL compute_domain CALL init_transfert CALL init_writefield + ! CALL allocate_field(sum_ne,field_T,type_real) ! CALL allocate_field_glo(sum_ne_glo,field_T,type_real) !