! ! $Id: etat0_netcdf.F 1322 2010-03-12 10:54:11Z emillour $ ! c c SUBROUTINE etat0_netcdf (interbar, masque) #ifdef CPP_EARTH USE startvar USE ioipsl USE dimphy USE control_mod USE infotrac USE fonte_neige_mod USE pbl_surface_mod USE phys_state_var_mod USE filtreg_mod use regr_lat_time_climoz_m, only: regr_lat_time_climoz use conf_phys_m, only: conf_phys #endif !#endif of #ifdef CPP_EARTH use netcdf, only: nf90_open, NF90_NOWRITE, nf90_close ! IMPLICIT NONE ! #include "dimensions.h" #include "paramet.h" ! ! ! INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, ! .KLON=KFDIA-KIDIA+1,KLEV=llm ! #ifdef CPP_EARTH #include "comgeom2.h" #include "comvert.h" #include "comconst.h" #include "indicesol.h" #include "dimsoil.h" #include "temps.h" #endif !#endif of #ifdef CPP_EARTH ! arguments: LOGICAL interbar REAL :: masque(iip1,jjp1) #ifdef CPP_EARTH ! local variables: REAL :: latfi(klon), lonfi(klon) REAL :: orog(iip1,jjp1), rugo(iip1,jjp1) REAL :: psol(iip1, jjp1), phis(iip1, jjp1) REAL :: p3d(iip1, jjp1, llm+1) REAL :: uvent(iip1, jjp1, llm) REAL :: vvent(iip1, jjm, llm) REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm) REAL :: qsat(iip1, jjp1, llm) REAL,ALLOCATABLE :: q3d(:, :, :,:) REAL :: tsol(klon), qsol(klon), sn(klon) !! REAL :: tsolsrf(klon,nbsrf) real qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) REAL :: albe(klon,nbsrf), evap(klon,nbsrf) REAL :: alblw(klon,nbsrf) REAL :: tsoil(klon,nsoilmx,nbsrf) REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf) REAL :: rugmer(klon) REAL :: qd(iip1, jjp1, llm) REAL :: run_off_lic_0(klon) ! declarations pour lecture glace de mer REAL :: rugv(klon) INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret INTEGER :: itaul(1), fid REAL :: lev(1), date REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic REAL :: flic_tmp(iip1, jjp1) REAL :: champint(iim, jjp1) ! CHARACTER(len=80) :: varname ! INTEGER :: i,j, ig, l, ji,ii1,ii2 REAL :: xpi ! REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1) REAL :: workvar(iip1,jjp1,llm) ! REAL :: prefkap, unskap ! real :: time_step,t_ops,t_wrt #include "comdissnew.h" #include "serre.h" #include "clesphys.h" INTEGER :: longcles PARAMETER ( longcles = 20 ) REAL :: clesphy0 ( longcles ) REAL :: p(iip1,jjp1,llm) INTEGER :: itau, iday REAL :: masse(iip1,jjp1,llm) REAL :: xpn,xps,xppn(iim),xpps(iim) real :: time REAL :: phi(ip1jmp1,llm) REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) REAL :: w(ip1jmp1,llm) REAL ::phystep CC REAL :: rugsrel(iip1*jjp1) REAL :: fder(klon) !! real zrel(iip1*jjp1),chmin,chmax !! CHARACTER(len=80) :: visu_file INTEGER :: visuid ! pour la lecture du fichier masque ocean integer :: nid_o2a logical :: couple = .false. INTEGER :: iml_omask, jml_omask REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp real, dimension(klon) :: ocemask_fi integer :: isst(klon-2) real zx_tmp_2d(iim,jjp1) REAL :: dummy logical :: ok_newmicro integer :: iflag_radia logical :: ok_journe, ok_mensuel, ok_instan, ok_hf logical :: ok_LES LOGICAL :: ok_ade, ok_aie, aerosol_couple, new_aod INTEGER :: flag_aerosol REAL :: bl95_b0, bl95_b1 real :: fact_cldcon, facttemps,ratqsbas,ratqshaut real :: tau_ratqs integer :: iflag_cldcon integer :: iflag_ratqs integer :: iflag_coupl integer :: iflag_clos integer :: iflag_wake integer :: iflag_thermals,nsplit_thermals real :: tau_thermals integer :: iflag_thermals_ed,iflag_thermals_optflux REAL :: solarlong0 real :: seuil_inversion real :: alp_offset integer read_climoz ! read ozone climatology C Allowed values are 0, 1 and 2 C 0: do not read an ozone climatology C 1: read a single ozone climatology that will be used day and night C 2: read two ozone climatologies, the average day and night C climatology and the daylight climatology ! ! Constantes ! pi = 4. * ATAN(1.) rad = 6371229. omeg = 4.* ASIN(1.)/(24.*3600.) g = 9.8 daysec = 86400. kappa = 0.2857143 cpp = 1004.70885 ! preff = 101325. pa = 50000. unskap = 1./kappa ! jmp1 = jjm + 1 ! ! Construct a grid ! ! CALL defrun_new(99,.TRUE.,clesphy0) CALL conf_gcm( 99, .TRUE. , clesphy0 ) call conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & & solarlong0,seuil_inversion, & & fact_cldcon, facttemps,ok_newmicro,iflag_radia, & & iflag_cldcon, & & iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & & ok_ade, ok_aie, aerosol_couple, & & flag_aerosol, new_aod, & & bl95_b0, bl95_b1, & & iflag_thermals,nsplit_thermals,tau_thermals, & & iflag_thermals_ed,iflag_thermals_optflux, & & iflag_coupl,iflag_clos,iflag_wake, read_climoz, & & alp_offset) ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value) co2_ppm0 = co2_ppm dtvr = daysec/REAL(day_step) print*,'dtvr',dtvr CALL iniconst() CALL inigeom() ! Initialisation pour traceurs call infotrac_init ALLOCATE(q3d(iip1, jjp1, llm, nqtot)) CALL inifilr() CALL phys_state_var_init(read_climoz) ! latfi(1) = ASIN(1.0) DO j = 2, jjm DO i = 1, iim latfi((j-2)*iim+1+i)= rlatu(j) ENDDO ENDDO latfi(klon) = - ASIN(1.0) ! lonfi(1) = 0.0 DO j = 2, jjm DO i = 1, iim lonfi((j-2)*iim+1+i) = rlonv(i) ENDDO ENDDO lonfi(klon) = 0.0 ! xpi = 2.0 * ASIN(1.0) DO ig = 1, klon latfi(ig) = latfi(ig) * 180.0 / xpi lonfi(ig) = lonfi(ig) * 180.0 / xpi ENDDO ! rlat(1) = ASIN(1.0) DO j = 2, jjm DO i = 1, iim rlat((j-2)*iim+1+i)= rlatu(j) ENDDO ENDDO rlat(klon) = - ASIN(1.0) ! rlon(1) = 0.0 DO j = 2, jjm DO i = 1, iim rlon((j-2)*iim+1+i) = rlonv(i) ENDDO ENDDO rlon(klon) = 0.0 ! xpi = 2.0 * ASIN(1.0) DO ig = 1, klon rlat(ig) = rlat(ig) * 180.0 / xpi rlon(ig) = rlon(ig) * 180.0 / xpi ENDDO ! C C En cas de simulation couplee, lecture du masque ocean issu du modele ocean C utilise pour calculer les poids et pour assurer l'adequation entre les C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque C a partir du fichier relief C write(*,*)'Essai de lecture masque ocean' iret = nf90_open("o2a.nc", NF90_NOWRITE, nid_o2a) if (iret .ne. 0) then write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve' write(*,*)'Run force' varname = 'masque' masque(:,:) = 0.0 CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, masque, $ 0.0, jjm ,rlonu,rlatv , interbar ) WRITE(*,*) 'MASQUE construit : Masque' WRITE(*,'(97I1)') nINT(masque(:,:)) call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq) WHERE (zmasq(1 : klon) .LT. EPSFRA) zmasq(1 : klon) = 0. END WHERE WHERE (1. - zmasq(1 : klon) .LT. EPSFRA) zmasq(1 : klon) = 1. END WHERE else couple = .true. iret = nf90_close(nid_o2a) call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp $ , nid_o2a) if (iml_omask /= iim .or. jml_omask /= jjp1) then write(*,*)'Dimensions non compatibles pour masque ocean' write(*,*)'iim = ',iim,' iml_omask = ',iml_omask write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask stop endif ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret) ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret) ALLOCATE(dlon_omask(iml_omask), stat=iret) ALLOCATE(dlat_omask(jml_omask), stat=iret) ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret) ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret) CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp $ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid) CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, $ ttm_tmp, 1, 1, ocetmp) CALL flinclo(fid) dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1) dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask) ocemask = ocetmp if (dlat_omask(1) < dlat_omask(jml_omask)) then do j = 1, jml_omask ocemask(:,j) = ocetmp(:,jml_omask-j+1) enddo endif C C passage masque ocean a la grille physique C write(*,*)'ocemask ' write(*,'(96i1)')int(ocemask) ocemask_fi(1) = ocemask(1,1) do j = 2, jjm do i = 1, iim ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j) enddo enddo ocemask_fi(klon) = ocemask(1,jjp1) zmasq = 1. - ocemask_fi endif call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) varname = 'relief' ! This line needs to be replaced by a call to restget to get the values in the restart file orog(:,:) = 0.0 CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, orog, $ 0.0 , jjm ,rlonu,rlatv , interbar, masque ) ! WRITE(*,*) 'OUT OF GET VARIABLE : Relief' ! WRITE(*,'(49I1)') INT(orog(:,:)) ! varname = 'rugosite' ! This line needs to be replaced by a call to restget to get the values in the restart file rugo(:,:) = 0.0 CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, rugo, $ 0.0 , jjm, rlonu,rlatv , interbar ) ! WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' ! WRITE(*,'(49I1)') INT(rugo(:,:)*10) ! C C on initialise les sous surfaces C pctsrf=0. c varname = 'psol' psol(:,:) = 0.0 CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, psol, $ 0.0 , jjm ,rlonu,rlatv , interbar ) ! ! Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM ! anyway. ! ! WRITE(*,*) 'PSOL :', psol(10,20) ! WRITE(*,*) ap(:), bp(:) CALL pression(ip1jmp1, ap, bp, psol, p3d) ! WRITE(*,*) 'P3D :', p3d(10,20,:) CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar) ! WRITE(*,*) 'PK:', pk(10,20,:) ! ! ! prefkap = preff ** kappa ! WRITE(*,*) 'unskap, cpp, preff :', unskap, cpp, preff DO l = 1, llm DO j=1,jjp1 DO i =1, iip1 pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap ENDDO ENDDO ENDDO ! ! WRITE(*,*) 'PLS :', pls(10,20,:) ! varname = 'surfgeo' phis(:,:) = 0.0 CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, phis, $ 0.0 , jjm ,rlonu,rlatv, interbar ) ! varname = 'u' uvent(:,:,:) = 0.0 CALL startget_dyn(varname, rlonu, rlatu, pls, workvar, uvent, 0., $ rlonv, rlatv, interbar ) ! varname = 'v' vvent(:,:,:) = 0.0 CALL startget_dyn(varname, rlonv, rlatv, pls(:, :jjm, :), . workvar(:, :jjm, :), vvent, 0., rlonu, rlatu(:jjm), interbar ) ! varname = 't' t3d(:,:,:) = 0.0 CALL startget_dyn(varname, rlonv, rlatu, pls, workvar, t3d, 0., $ rlonu, rlatv , interbar ) ! WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)), . maxval(t3d(:,:,:)) varname = 'tpot' tpot(:,:,:) = 0.0 CALL startget_dyn(varname, rlonv, rlatu, pls, pk, tpot, 0., rlonu, $ rlatv, interbar) ! WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)), . maxval(t3d(:,:,:)) WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)), . maxval(pls(:,:,:)) c Calcul de l'humidite a saturation print*,'avant q_sat' call q_sat(llm*jjp1*iip1,t3d,pls,qsat) print*,'apres q_sat' WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)), . maxval(qsat(:,:,:)) ! CC WRITE(*,*) 'QSAT :', qsat(10,20,:) ! varname = 'q' qd(:,:,:) = 0.0 q3d(:,:,:,:) = 0.0 WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)), . maxval(qsat(:,:,:)) CALL startget_dyn(varname, rlonv, rlatu, pls, qsat, qd, 0., rlonu, $ rlatv , interbar ) q3d(:,:,:,1) = qd(:,:,:) ! ! Ozone climatology: if (read_climoz >= 1) call regr_lat_time_climoz(read_climoz) varname = 'tsol' ! This line needs to be replaced by a call to restget to get the values in the restart file tsol(:) = 0.0 CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon, $ tsol, 0.0, jjm, rlonu, rlatv , interbar ) ! WRITE(*,*) 'TSOL construit :' ! WRITE(*,'(48I3)') INT(TSOL(2:klon)-273) ! varname = 'qsol' qsol(:) = 0.0 CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon, $ qsol, 0.0, jjm, rlonu, rlatv , interbar ) ! varname = 'snow' sn(:) = 0.0 CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon, sn, $ 0.0, jjm, rlonu, rlatv , interbar ) ! varname = 'rads' radsol(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,radsol, $ 0.0, jjm, rlonu, rlatv , interbar ) ! varname = 'rugmer' rugmer(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer, $ 0.0, jjm, rlonu, rlatv , interbar ) ! ! varname = 'agesno' ! agesno(:) = 0.0 ! CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0, ! . jjm, rlonu, rlatv , interbar ) varname = 'zmea' zmea(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0, . jjm, rlonu, rlatv , interbar ) varname = 'zstd' zstd(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0, . jjm, rlonu, rlatv , interbar ) varname = 'zsig' zsig(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0, . jjm, rlonu, rlatv , interbar ) varname = 'zgam' zgam(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0, . jjm, rlonu, rlatv , interbar ) varname = 'zthe' zthe(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0, . jjm, rlonu, rlatv , interbar ) varname = 'zpic' zpic(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0, . jjm, rlonu, rlatv , interbar ) varname = 'zval' zval(:) = 0.0 CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0, . jjm, rlonu, rlatv , interbar ) c cc rugsrel(:) = 0.0 cc IF(ok_orodr) THEN cc DO i = 1, iip1* jjp1 cc rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. ) cc ENDDO cc ENDIF C C lecture du fichier glace de terre pour fixer la fraction de terre C et de glace de terre C CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp $ , fid) ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret) ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret) ALLOCATE(dlon_lic(iml_lic), stat=iret) ALLOCATE(dlat_lic(jml_lic), stat=iret) ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret) CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp $ , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp $ , 1, 1, fraclic) CALL flinclo(fid) C C interpolation sur la grille T du modele C WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ', $ iml_lic, jml_lic c C sil les coordonnees sont en degres, on les transforme C IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) ) THEN lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180. ENDIF IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180. ENDIF dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1) dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic) C CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic $ ,iim, jjp1, $ rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1)) cx$$$ flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1) flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1) C C passage sur la grille physique C CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, $ pctsrf(1:klon, is_lic)) C adequation avec le maque terre/mer c zmasq(157) = 0. WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA ) pctsrf(1 : klon, is_lic) = 0. END WHERE WHERE (zmasq( 1 : klon) .LT. EPSFRA) pctsrf(1 : klon, is_lic) = 0. END WHERE pctsrf(1 : klon, is_ter) = zmasq(1 : klon) DO ji = 1, klon IF (zmasq(ji) .GT. EPSFRA) THEN IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN pctsrf(ji, is_lic) = zmasq(ji) pctsrf(ji, is_ter) = 0. ELSE pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic) IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN pctsrf(ji,is_ter) = 0. pctsrf(ji, is_lic) = zmasq(ji) ENDIF ENDIF ENDIF END DO C C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0) C pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon)) WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA) pctsrf(1 : klon, is_oce) = 0. END WHERE if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon) isst = 0 where (pctsrf(2:klon-1,is_oce) >0.) isst = 1 C C verif que somme des sous surface = 1 C ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf),dim=2))-1.0) $ .GT. EPSFRA) IF (ji .NE. 0) THEN WRITE(*,*) 'pb repartition sous maille pour ',ji,' points' ENDIF ! where (pctsrf(1:klon, is_ter) >= .5) ! pctsrf(1:klon, is_ter) = 1. ! pctsrf(1:klon, is_oce) = 0. ! pctsrf(1:klon, is_sic) = 0. ! pctsrf(1:klon, is_lic) = 0. ! zmasq = 1. ! endwhere ! where (pctsrf(1:klon, is_lic) >= .5) ! pctsrf(1:klon, is_ter) = 0. ! pctsrf(1:klon, is_oce) = 0. ! pctsrf(1:klon, is_sic) = 0. ! pctsrf(1:klon, is_lic) = 1. ! zmasq = 1. ! endwhere ! where (pctsrf(1:klon, is_oce) >= .5) ! pctsrf(1:klon, is_ter) = 0. ! pctsrf(1:klon, is_oce) = 1. ! pctsrf(1:klon, is_sic) = 0. ! pctsrf(1:klon, is_lic) = 0. ! zmasq = 0. ! endwhere ! where (pctsrf(1:klon, is_sic) >= .5) ! pctsrf(1:klon, is_ter) = 0. ! pctsrf(1:klon, is_oce) = 0. ! pctsrf(1:klon, is_sic) = 1. ! pctsrf(1:klon, is_lic) = 0. ! zmasq = 0. ! endwhere ! call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) C C verif que somme des sous surface = 1 C ! ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 ) ! $ .GT. EPSFRA) ! IF (ji .NE. 0) THEN ! WRITE(*,*) 'pb repartition sous maille pour ',ji,' points' ! ENDIF CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d) write(*,*)'zmasq = ' write(*,'(96i1)')nint(zx_tmp_2d) call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) WRITE(*,*) 'MASQUE construit : Masque' WRITE(*,'(97I1)') nINT(masque(:,:)) C Calcul intermediaire c CALL massdair( p3d, masse ) c print *,' ALPHAX ',alphax DO l = 1, llm DO i = 1, iim xppn(i) = aire( i, 1 ) * masse( i , 1 , l ) xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l ) ENDDO xpn = SUM(xppn)/apoln xps = SUM(xpps)/apols DO i = 1, iip1 masse( i , 1 , l ) = xpn masse( i , jjp1 , l ) = xps ENDDO ENDDO q3d(iip1,:,:,:) = q3d(1,:,:,:) phis(iip1,:) = phis(1,:) C Ecriture CALL inidissip( lstardis, nitergdiv, nitergrot, niterh , * tetagdiv, tetagrot , tetatemp ) print*,'sortie inidissip' itau = 0 itau_dyn = 0 itau_phy = 0 iday = dayref +itau/day_step time = real(itau-(iday-dayref)*day_step)/day_step c IF(time.GT.1) THEN time = time - 1 iday = iday + 1 ENDIF day_ref = dayref annee_ref = anneeref CALL geopot ( ip1jmp1, tpot , pk , pks, phis , phi ) print*,'sortie geopot' CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis , * phi,w, pbaru,pbarv,time+iday-dayref ) print*,'sortie caldyn0' CALL dynredem0("start.nc",dayref,phis) print*,'sortie dynredem0' CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,masse , . psol) print*,'sortie dynredem1' C C Ecriture etat initial physique C write(*,*)'phystep ',dtvr,iphysiq,nbapp_rad phystep = dtvr * REAL(iphysiq) radpas = NINT (86400./phystep/ REAL(nbapp_rad) ) write(*,*)'phystep =', phystep, radpas cIM : lecture de co2_ppm & solaire ds physiq.def c co2_ppm = 348.0 c solaire = 1365.0 c c Initialisation c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs c ftsol(:,is_ter) = tsol ftsol(:,is_lic) = tsol ftsol(:,is_oce) = tsol ftsol(:,is_sic) = tsol snsrf(:,is_ter) = sn snsrf(:,is_lic) = sn snsrf(:,is_oce) = sn snsrf(:,is_sic) = sn falb1(:,is_ter) = 0.08 falb1(:,is_lic) = 0.6 falb1(:,is_oce) = 0.5 falb1(:,is_sic) = 0.6 falb2 = falb1 evap(:,:) = 0. qsolsrf(:,is_ter) = 150 qsolsrf(:,is_lic) = 150 qsolsrf(:,is_oce) = 150. qsolsrf(:,is_sic) = 150. do i = 1, nbsrf do j = 1, nsoilmx tsoil(:,j,i) = tsol enddo enddo rain_fall = 0.; snow_fall = 0. solsw = 165. sollw = -53. t_ancien = 273.15 q_ancien = 0. agesno = 0. c frugs(1:klon,is_oce) = rugmer(1:klon) frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0) frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0) frugs(1:klon,is_sic) = 0.001 fder = 0.0 clwcon = 0.0 rnebcon = 0.0 ratqs = 0.0 run_off_lic_0 = 0.0 rugoro = 0.0 c c Avant l'appel a phyredem, on initialize les modules de surface c avec les valeurs qui vont etre ecrit dans startphy.nc c dummy = 1.0 pbl_tke(:,:,:) = 1.e-8 zmax0(:) = 40. f0(:) = 1.e-5 ema_work1(:,:) = 0. ema_work2(:,:) = 0. wake_deltat(:,:) = 0. wake_deltaq(:,:) = 0. wake_s(:) = 0. wake_cstar(:) = 0. wake_fip(:) = 0. call fonte_neige_init(run_off_lic_0) call pbl_surface_init(qsol, fder, snsrf, qsolsrf, $ evap, frugs, agesno, tsoil) call phyredem("startphy.nc") C Sortie Visu pour les champs dynamiques cc if (1.eq.0 ) then cc print*,'sortie visu' cc time_step = 1. cc t_ops = 2. cc t_wrt = 2. cc itau = 2. cc visu_file='Etat0_visu.nc' cc CALL initdynav(visu_file,dayref,anneeref,time_step, cc . t_ops, t_wrt, visuid) cc CALL writedynav(visuid, itau,vvent , cc . uvent,tpot,pk,phi,q3d,masse,psol,phis) cc else print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0' cc endif print*,'entree histclo' CALL histclo #endif !#endif of #ifdef CPP_EARTH RETURN ! END SUBROUTINE etat0_netcdf