SUBROUTINE etat0_netcdf USE startvar USE ioipsl ! IMPLICIT NONE ! #include "netcdf.inc" #include "dimensions.h" #include "paramet.h" ! ! c INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, c .KLON=KFDIA-KIDIA+1,KLEV=llm ! #include "comgeom2.h" #include "comvert.h" #include "comconst.h" #include "indicesol.h" #include "dimphy.h" #include "dimsoil.h" ! REAL :: latfi(klon), lonfi(klon) REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1), . 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 :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm) REAL :: tsol(klon), qsol(klon), sn(klon) REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) REAL :: albe(klon,nbsrf), evap(klon,nbsrf) REAL :: tsoil(klon,nsoilmx,nbsrf) REAL :: radsol(klon),rain_fall(klon), snow_fall(klon) REAL :: solsw(klon), sollw(klon) REAL :: deltat(klon), frugs(klon,nbsrf), agesno(klon),rugmer(klon) REAL :: zmea(iip1*jjp1), zstd(iip1*jjp1) REAL :: zsig(iip1*jjp1), zgam(iip1*jjp1), zthe(iip1*jjp1) REAL :: zpic(iip1*jjp1), zval(iip1*jjp1), rugsrel(iip1*jjp1) REAL :: qd(iip1, jjp1, llm) REAL :: pctsrf(klon, nbsrf) REAL :: t_ancien(klon,klev), q_ancien(klon,klev) ! ! 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, dt 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*80 :: varname ! INTEGER :: i,j, ig, l, ji 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 :: q_sat EXTERNAL q_sat real :: time_step,t_ops,t_wrt #include "comdissnew.h" #include "control.h" #include "serre.h" #include "clesph0.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,co2_ppm,solaire INTEGER :: radpas CHARACTER*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) ! ! 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. unskap = 1./kappa ! jmp1 = jjm + 1 ! ! Construct a grid ! ! CALL defrun_new(99,.TRUE.,clesphy0) CALL conf_gcm( 99, .TRUE. , clesphy0 ) dtvr = daysec/FLOAT(day_step) print*,'dtvr',dtvr CALL inicons0() CALL inigeom() ! CALL inifilr() ! 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 ! 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(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0) ! 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(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0) ! WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' WRITE(*,'(49I1)') INT(rugo(:,:)*10) ! varname = 'masque' ! This line needs to be replaced by a call to restget to get the values in the restart file masque(:,:) = 0.0 CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0) ! WRITE(*,*) 'MASQUE construit : Masque' WRITE(*,'(49I1)') INT(masque(:,:)) ! ! C C on initialise les sous surfaces C pctsrf=0. !cree le masque a partir du fichier relief varname = 'zmasq' zmasq(:) = 0. CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0) 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 WRITE(*,*)zmasq varname = 'psol' psol(:,:) = 0.0 CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0) ! ! 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(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0) write(*,*) 'Phis = ' write(*,*)phis ! varname = 'u' uvent(:,:,:) = 0.0 CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls, . workvar, uvent, 0.0) ! varname = 'v' vvent(:,:,:) = 0.0 CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls, . workvar, vvent, 0.0) ! varname = 't' t3d(:,:,:) = 0.0 CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, . workvar, t3d, 0.0) ! WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)), . maxval(t3d(:,:,:)) varname = 'tpot' tpot(:,:,:) = 0.0 CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, . pk, tpot, 0.0) ! WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)), . maxval(t3d(:,:,:)) WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)), . maxval(pls(:,:,:)) DO l = 1, llm DO j=1,jjp1 DO i =1, iip1-1 qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. ) ENDDO qsat(iip1,j,l) = qsat(1,j,l) ENDDO ENDDO WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)), . maxval(qsat(:,:,:)) ! WRITE(*,*) 'QSAT :', qsat(10,20,:) ! varname = 'q' q3d(:,:,:,:) = 0.0 qd(:,:,:) = 0.0 q3d(:,:,:,:) = 0.0 WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)), . maxval(qsat(:,:,:)) CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, . qsat, qd, 0.0) q3d(:,:,:,1) = qd(:,:,:) ! 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(varname, iip1, jjp1, rlonv, rlatu, klon, tsol,0.0) ! WRITE(*,*) 'TSOL construit :' WRITE(*,'(48I3)') INT(TSOL(2:klon)-273) ! varname = 'qsol' qsol(:) = 0.0 CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol,0.0) ! varname = 'snow' sn(:) = 0.0 CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn,0.0) ! varname = 'rads' radsol(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0) ! varname = 'deltat' deltat(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0) ! varname = 'rugmer' rugmer(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0) ! varname = 'agsno' agesno(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0) varname = 'zmea' zmea(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0) varname = 'zstd' zstd(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0) varname = 'zsig' zsig(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0) varname = 'zgam' zgam(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0) varname = 'zthe' zthe(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0) varname = 'zpic' zpic(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0) varname = 'zval' zval(:) = 0.0 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0) rugsrel(:) = 0.0 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 C write(*,*)'Essai de lecture masque ocean' iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a) if (iret .ne. 0) then write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve' write(*,*)'Run force' else couple = .true. iret = nf_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 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 isst = 0 where (ocemask_fi(2:klon-1) >0.) isst = 1 write(45,'(72i1)')isst 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)) c$$$ 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 write(46,'(72i1)')isst 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 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 iday = dayref +itau/day_step time = FLOAT(itau-(iday-dayref)*day_step)/day_step c IF(time.GT.1) THEN time = time - 1 iday = iday + 1 ENDIF 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,anneeref,phis,nqmx) print*,'sortie dynredem0' CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse , . psol) print*,'sortie dynredem1' C C Ecriture etat initial physique C phystep = dtvr * FLOAT(iphysiq) radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) ) co2_ppm = 330.0 solaire = 1370.0 c call physdem(lonfi, latfi, phystep,radpas,co2_ppm, c . solaire,tsol, qsol, c . sn, radsol, deltat, rugmer, c . agesno, zmea, zstd, zsig, c . zgam, zthe, zpic, zval, c . rugsrel) c c Initialisation c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs c tsolsrf(:,is_ter) = tsol tsolsrf(:,is_lic) = tsol tsolsrf(:,is_oce) = tsol tsolsrf(:,is_sic) = tsol snsrf(:,is_ter) = sn snsrf(:,is_lic) = sn snsrf(:,is_oce) = sn snsrf(:,is_sic) = sn albe(:,is_ter) = 0.08 albe(:,is_lic) = 0.6 albe(:,is_oce) = 0.5 albe(:,is_sic) = 0.6 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. deltat = 0. 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 call physdem("startphy.nc",phystep,radpas, co2_ppm, solaire, $ latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf, $ albe, evap, rain_fall, snow_fall, solsw, sollw, $ radsol, frugs, agesno, $ zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, $ t_ancien, q_ancien) C Sortie Visu pour les champs dynamiques print*,'sortie visu' time_step = 1. t_ops = 2. t_wrt = 2. itau = 2. visu_file='Etat0_visu.nc' CALL initdynav(visu_file,dayref,anneeref,time_step, . t_ops, t_wrt, nqmx, visuid) CALL writedynav(visuid, nqmx, itau,vvent , . uvent,tpot,pk,phi,q3d,masse,psol,phis) print*,'entree histclo' CALL histclo RETURN ! END SUBROUTINE etat0_netcdf