SUBROUTINE etat0_netcdf USE startvar USE ioipsl ! IMPLICIT NONE ! #include "dimensions.h" #include "paramet.h" ! ! INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, .KLON=KFDIA-KIDIA+1,KLEV=llm ! #include "comgeom2.h" #include "comvert.h" #include "comconst.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), radsol(klon) REAL :: deltat(klon), rugmer(klon), agesno(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) ! CHARACTER*80 :: varname ! INTEGER :: i,j, ig, l 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 ! ! 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) 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(:,:)) ! ! 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) ! 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' qd(:,:,:) = 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 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 call physdem(lonfi, latfi, phystep,radpas,co2_ppm, . solaire,tsol, qsol, . sn, radsol, deltat, rugmer, . agesno, zmea, zstd, zsig, . zgam, zthe, zpic, zval, . rugsrel) 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