Changeset 1319 for LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90
- Timestamp:
- Feb 23, 2010, 10:29:54 PM (14 years ago)
- File:
-
- 1 moved
-
LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90 (moved) (moved from LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90
r1318 r1319 1 1 ! 2 ! $ Header$2 ! $Id$ 3 3 ! 4 c 5 c 6 SUBROUTINE etat0_netcdf (interbar, masque) 7 #ifdef CPP_EARTH 8 USE startvar 9 USE ioipsl 10 USE dimphy 11 USE infotrac 12 USE fonte_neige_mod 13 USE pbl_surface_mod 14 USE phys_state_var_mod 15 USE filtreg_mod 16 use regr_lat_time_climoz_m, only: regr_lat_time_climoz 17 use conf_phys_m, only: conf_phys 4 !------------------------------------------------------------------------------- 5 ! 6 SUBROUTINE etat0_netcdf(ib, masque, letat0) 7 ! 8 !------------------------------------------------------------------------------- 9 ! Purpose: Creates initial states 10 !------------------------------------------------------------------------------- 11 ! Note: This routine is designed to work for Earth 12 !------------------------------------------------------------------------------- 13 #ifdef CPP_EARTH 14 USE startvar 15 USE ioipsl 16 USE dimphy 17 USE infotrac 18 USE fonte_neige_mod 19 USE pbl_surface_mod 20 USE phys_state_var_mod 21 USE filtreg_mod 22 USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz 23 USE conf_phys_m, ONLY: conf_phys 24 USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR 18 25 #endif 19 !#endif of #ifdef CPP_EARTH 20 use netcdf, only: nf90_open, NF90_NOWRITE, nf90_close 21 ! 22 IMPLICIT NONE 23 ! 26 IMPLICIT NONE 27 !------------------------------------------------------------------------------- 28 ! Arguments: 24 29 #include "dimensions.h" 25 30 #include "paramet.h" 26 ! 27 ! 28 ! INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, 29 ! .KLON=KFDIA-KIDIA+1,KLEV=llm 30 ! 31 #ifdef CPP_EARTH 31 #include "iniprint.h" 32 LOGICAL, INTENT(IN) :: ib ! barycentric interpolat. 33 REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask 34 LOGICAL, INTENT(IN) :: letat0 ! F: masque only required 35 #ifndef CPP_EARTH 36 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' 37 #else 38 !------------------------------------------------------------------------------- 39 ! Local variables: 32 40 #include "comgeom2.h" 33 41 #include "comvert.h" … … 36 44 #include "dimsoil.h" 37 45 #include "temps.h" 38 #endif 39 !#endif of #ifdef CPP_EARTH 40 ! arguments: 41 LOGICAL interbar 42 REAL :: masque(iip1,jjp1) 43 44 #ifdef CPP_EARTH 45 ! local variables: 46 REAL :: latfi(klon), lonfi(klon) 47 REAL :: orog(iip1,jjp1), rugo(iip1,jjp1) 48 REAL :: psol(iip1, jjp1), phis(iip1, jjp1) 49 REAL :: p3d(iip1, jjp1, llm+1) 50 REAL :: uvent(iip1, jjp1, llm) 51 REAL :: vvent(iip1, jjm, llm) 52 REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm) 53 REAL :: qsat(iip1, jjp1, llm) 54 REAL,ALLOCATABLE :: q3d(:, :, :,:) 55 REAL :: tsol(klon), qsol(klon), sn(klon) 56 !! REAL :: tsolsrf(klon,nbsrf) 57 real qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 58 REAL :: albe(klon,nbsrf), evap(klon,nbsrf) 59 REAL :: alblw(klon,nbsrf) 60 REAL :: tsoil(klon,nsoilmx,nbsrf) 61 REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf) 62 REAL :: rugmer(klon) 63 REAL :: qd(iip1, jjp1, llm) 64 REAL :: run_off_lic_0(klon) 65 ! declarations pour lecture glace de mer 66 REAL :: rugv(klon) 67 INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret 68 INTEGER :: itaul(1), fid 69 REAL :: lev(1), date 70 REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic 71 REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic 72 REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic 73 REAL :: flic_tmp(iip1, jjp1) 74 REAL :: champint(iim, jjp1) 75 ! 76 77 CHARACTER(len=80) :: varname 78 ! 79 INTEGER :: i,j, ig, l, ji,ii1,ii2 80 REAL :: xpi 81 ! 82 REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) 83 REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1) 84 REAL :: workvar(iip1,jjp1,llm) 85 ! 86 REAL :: prefkap, unskap 87 ! 88 real :: time_step,t_ops,t_wrt 46 REAL, DIMENSION(klon) :: tsol, qsol 47 REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0 48 REAL, DIMENSION(iip1,jjp1) :: orog, rugo, psol, phis 49 REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d 50 REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd 51 REAL, DIMENSION(iip1,jjm ,llm) :: vvent 52 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: q3d 53 REAL, DIMENSION(klon,nbsrf) :: qsolsrf, snsrf, evap 54 REAL, DIMENSION(klon,nbsrf) :: frugs, agesno 55 REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil 56 57 !--- Local variables for sea-ice reading: 58 INTEGER :: iml_lic, jml_lic, llm_tmp 59 INTEGER :: ttm_tmp, iret, fid 60 INTEGER, DIMENSION(1) :: itaul 61 REAL, DIMENSION(1) :: lev 62 REAL :: date 63 REAL, DIMENSION(:,:), ALLOCATABLE :: lon_lic, lat_lic, fraclic 64 REAL, DIMENSION(:), ALLOCATABLE :: dlon_lic, dlat_lic 65 REAL, DIMENSION(iip1,jjp1) :: flic_tmp 66 67 !--- Misc 68 CHARACTER(LEN=80) :: x, fmt 69 INTEGER :: i, j, l, ji 70 REAL, DIMENSION(iip1,jjp1,llm) :: alpha, beta, pk, pls, y 71 REAL, DIMENSION(ip1jmp1) :: pks 89 72 90 73 #include "comdissnew.h" … … 93 76 #include "clesphys.h" 94 77 95 INTEGER :: longcles 96 PARAMETER ( longcles = 20 ) 97 REAL :: clesphy0 ( longcles ) 98 REAL :: p(iip1,jjp1,llm) 99 INTEGER :: itau, iday 100 REAL :: masse(iip1,jjp1,llm) 101 REAL :: xpn,xps,xppn(iim),xpps(iim) 102 real :: time 103 REAL :: phi(ip1jmp1,llm) 104 REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 105 REAL :: w(ip1jmp1,llm) 106 REAL ::phystep 107 CC REAL :: rugsrel(iip1*jjp1) 108 REAL :: fder(klon) 109 !! real zrel(iip1*jjp1),chmin,chmax 110 111 !! CHARACTER(len=80) :: visu_file 112 INTEGER :: visuid 113 114 ! pour la lecture du fichier masque ocean 115 integer :: nid_o2a 116 logical :: couple = .false. 117 INTEGER :: iml_omask, jml_omask 118 REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask 119 REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask 120 REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp 121 real, dimension(klon) :: ocemask_fi 122 integer :: isst(klon-2) 123 real zx_tmp_2d(iim,jjp1) 124 125 REAL :: dummy 126 127 logical :: ok_newmicro 128 integer :: iflag_radia 129 logical :: ok_journe, ok_mensuel, ok_instan, ok_hf 130 logical :: ok_LES 131 LOGICAL :: ok_ade, ok_aie, aerosol_couple, new_aod 132 INTEGER :: flag_aerosol 133 REAL :: bl95_b0, bl95_b1 134 real :: fact_cldcon, facttemps,ratqsbas,ratqshaut 135 real :: tau_ratqs 136 integer :: iflag_cldcon 137 integer :: iflag_ratqs 138 integer :: iflag_coupl 139 integer :: iflag_clos 140 integer :: iflag_wake 141 integer :: iflag_thermals,nsplit_thermals 142 real :: tau_thermals 143 integer :: iflag_thermals_ed,iflag_thermals_optflux 144 REAL :: solarlong0 145 real :: seuil_inversion 146 147 integer read_climoz ! read ozone climatology 148 C Allowed values are 0, 1 and 2 149 C 0: do not read an ozone climatology 150 C 1: read a single ozone climatology that will be used day and night 151 C 2: read two ozone climatologies, the average day and night 152 C climatology and the daylight climatology 153 154 ! 155 ! Constantes 156 ! 157 pi = 4. * ATAN(1.) 158 rad = 6371229. 159 omeg = 4.* ASIN(1.)/(24.*3600.) 160 g = 9.8 161 daysec = 86400. 162 kappa = 0.2857143 163 cpp = 1004.70885 164 ! 165 preff = 101325. 166 pa = 50000. 167 unskap = 1./kappa 168 ! 169 jmp1 = jjm + 1 170 ! 171 ! Construct a grid 172 ! 173 174 ! CALL defrun_new(99,.TRUE.,clesphy0) 175 CALL conf_gcm( 99, .TRUE. , clesphy0 ) 176 call conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & 177 & solarlong0,seuil_inversion, & 178 & fact_cldcon, facttemps,ok_newmicro,iflag_radia, & 179 & iflag_cldcon, & 180 & iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & 181 & ok_ade, ok_aie, aerosol_couple, & 182 & flag_aerosol, new_aod, & 183 & bl95_b0, bl95_b1, & 184 & iflag_thermals,nsplit_thermals,tau_thermals, & 185 & iflag_thermals_ed,iflag_thermals_optflux, & 186 & iflag_coupl,iflag_clos,iflag_wake, read_climoz ) 78 REAL, DIMENSION(iip1,jjp1,llm) :: masse 79 INTEGER :: itau, iday 80 REAL :: xpn, xps, time, phystep 81 REAL, DIMENSION(iim) :: xppn, xpps 82 REAL, DIMENSION(ip1jmp1,llm) :: pbaru, phi, w 83 REAL, DIMENSION(ip1jm ,llm) :: pbarv 84 REAL, DIMENSION(klon) :: fder 85 86 !--- Local variables for ocean mask reading: 87 INTEGER :: nid_o2a, iml_omask, jml_omask 88 LOGICAL :: couple=.FALSE. 89 REAL, DIMENSION(:,:), ALLOCATABLE :: lon_omask, lat_omask, ocemask, ocetmp 90 REAL, DIMENSION(:), ALLOCATABLE :: dlon_omask,dlat_omask 91 REAL, DIMENSION(klon) :: ocemask_fi 92 INTEGER, DIMENSION(klon-2) :: isst 93 REAL, DIMENSION(iim,jjp1) :: zx_tmp_2d 94 REAL :: dummy 95 LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf 96 LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod 97 INTEGER :: iflag_radia, flag_aerosol 98 REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut 99 REAL :: tau_ratqs 100 INTEGER :: iflag_cldcon, iflag_ratqs, iflag_coupl, iflag_clos, iflag_wake 101 INTEGER :: iflag_thermals, nsplit_thermals 102 INTEGER :: iflag_thermals_ed, iflag_thermals_optflux 103 REAL :: tau_thermals, solarlong0, seuil_inversion 104 INTEGER :: read_climoz ! read ozone climatology 105 ! Allowed values are 0, 1 and 2 106 ! 0: do not read an ozone climatology 107 ! 1: read a single ozone climatology that will be used day and night 108 ! 2: read two ozone climatologies, the average day and night 109 ! climatology and the daylight climatology 110 !------------------------------------------------------------------------------- 111 !--- Constants 112 pi = 4. * ATAN(1.) 113 rad = 6371229. 114 daysec = 86400. 115 omeg = 2.*pi/daysec 116 g = 9.8 117 kappa = 0.2857143 118 cpp = 1004.70885 119 preff = 101325. 120 pa = 50000. 121 jmp1 = jjm + 1 122 123 !--- CONSTRUCT A GRID 124 CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & 125 solarlong0,seuil_inversion, & 126 fact_cldcon, facttemps,ok_newmicro,iflag_radia, & 127 iflag_cldcon, & 128 iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & 129 ok_ade, ok_aie, aerosol_couple, & 130 flag_aerosol, new_aod, & 131 bl95_b0, bl95_b1, & 132 iflag_thermals,nsplit_thermals,tau_thermals, & 133 iflag_thermals_ed,iflag_thermals_optflux, & 134 iflag_coupl,iflag_clos,iflag_wake, read_climoz ) 187 135 188 136 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value) 189 co2_ppm0 = co2_ppm 190 191 dtvr = daysec/FLOAT(day_step) 192 print*,'dtvr',dtvr 193 194 CALL iniconst() 195 CALL inigeom() 196 197 ! Initialisation pour traceurs 198 call infotrac_init 199 ALLOCATE(q3d(iip1, jjp1, llm, nqtot)) 200 201 CALL inifilr() 202 CALL phys_state_var_init(read_climoz) 203 ! 204 latfi(1) = ASIN(1.0) 205 DO j = 2, jjm 206 DO i = 1, iim 207 latfi((j-2)*iim+1+i)= rlatu(j) 208 ENDDO 209 ENDDO 210 latfi(klon) = - ASIN(1.0) 211 ! 212 lonfi(1) = 0.0 213 DO j = 2, jjm 214 DO i = 1, iim 215 lonfi((j-2)*iim+1+i) = rlonv(i) 216 ENDDO 217 ENDDO 218 lonfi(klon) = 0.0 219 ! 220 xpi = 2.0 * ASIN(1.0) 221 DO ig = 1, klon 222 latfi(ig) = latfi(ig) * 180.0 / xpi 223 lonfi(ig) = lonfi(ig) * 180.0 / xpi 224 ENDDO 225 ! 226 rlat(1) = ASIN(1.0) 227 DO j = 2, jjm 228 DO i = 1, iim 229 rlat((j-2)*iim+1+i)= rlatu(j) 230 ENDDO 231 ENDDO 232 rlat(klon) = - ASIN(1.0) 233 ! 234 rlon(1) = 0.0 235 DO j = 2, jjm 236 DO i = 1, iim 237 rlon((j-2)*iim+1+i) = rlonv(i) 238 ENDDO 239 ENDDO 240 rlon(klon) = 0.0 241 ! 242 xpi = 2.0 * ASIN(1.0) 243 DO ig = 1, klon 244 rlat(ig) = rlat(ig) * 180.0 / xpi 245 rlon(ig) = rlon(ig) * 180.0 / xpi 246 ENDDO 247 ! 248 249 250 251 C 252 C En cas de simulation couplee, lecture du masque ocean issu du modele ocean 253 C utilise pour calculer les poids et pour assurer l'adequation entre les 254 C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque 255 C a partir du fichier relief 256 C 257 258 write(*,*)'Essai de lecture masque ocean' 259 iret = nf90_open("o2a.nc", NF90_NOWRITE, nid_o2a) 260 if (iret .ne. 0) then 261 write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve' 262 write(*,*)'Run force' 263 varname = 'masque' 264 masque(:,:) = 0.0 265 CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0, 266 , jjm ,rlonu,rlatv , interbar ) 267 WRITE(*,*) 'MASQUE construit : Masque' 268 WRITE(*,'(97I1)') nINT(masque(:,:)) 269 call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq) 270 WHERE (zmasq(1 : klon) .LT. EPSFRA) 271 zmasq(1 : klon) = 0. 272 END WHERE 273 WHERE (1. - zmasq(1 : klon) .LT. EPSFRA) 274 zmasq(1 : klon) = 1. 275 END WHERE 276 else 277 couple = .true. 278 iret = nf90_close(nid_o2a) 279 call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp 280 $ , nid_o2a) 281 if (iml_omask /= iim .or. jml_omask /= jjp1) then 282 write(*,*)'Dimensions non compatibles pour masque ocean' 283 write(*,*)'iim = ',iim,' iml_omask = ',iml_omask 284 write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask 285 stop 286 endif 287 ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret) 288 ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret) 289 ALLOCATE(dlon_omask(iml_omask), stat=iret) 290 ALLOCATE(dlat_omask(jml_omask), stat=iret) 291 ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret) 292 ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret) 293 CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp 294 $ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid) 295 CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, 296 $ ttm_tmp, 1, 1, ocetmp) 297 CALL flinclo(fid) 298 dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1) 299 dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask) 300 ocemask = ocetmp 301 if (dlat_omask(1) < dlat_omask(jml_omask)) then 302 do j = 1, jml_omask 303 ocemask(:,j) = ocetmp(:,jml_omask-j+1) 304 enddo 305 endif 306 C 307 C passage masque ocean a la grille physique 308 C 309 write(*,*)'ocemask ' 310 write(*,'(96i1)')int(ocemask) 311 ocemask_fi(1) = ocemask(1,1) 312 do j = 2, jjm 313 do i = 1, iim 314 ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j) 315 enddo 316 enddo 317 ocemask_fi(klon) = ocemask(1,jjp1) 318 zmasq = 1. - ocemask_fi 319 endif 320 321 call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) 322 323 varname = 'relief' 324 ! This line needs to be replaced by a call to restget to get the values in the restart file 325 orog(:,:) = 0.0 326 CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 , 327 , jjm ,rlonu,rlatv , interbar, masque ) 328 ! 329 WRITE(*,*) 'OUT OF GET VARIABLE : Relief' 330 ! WRITE(*,'(49I1)') INT(orog(:,:)) 331 ! 332 varname = 'rugosite' 333 ! This line needs to be replaced by a call to restget to get the values in the restart file 334 rugo(:,:) = 0.0 335 CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 , 336 , jjm, rlonu,rlatv , interbar ) 337 ! 338 WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' 339 ! WRITE(*,'(49I1)') INT(rugo(:,:)*10) 340 ! 341 C 342 C on initialise les sous surfaces 343 C 344 pctsrf=0. 345 c 346 varname = 'psol' 347 psol(:,:) = 0.0 348 CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 , 349 , jjm ,rlonu,rlatv , interbar ) 350 ! 351 ! Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM 352 ! anyway. 353 ! 354 ! WRITE(*,*) 'PSOL :', psol(10,20) 355 ! WRITE(*,*) ap(:), bp(:) 356 CALL pression(ip1jmp1, ap, bp, psol, p3d) 357 ! WRITE(*,*) 'P3D :', p3d(10,20,:) 358 CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar) 359 ! WRITE(*,*) 'PK:', pk(10,20,:) 360 ! 361 ! 362 ! 363 prefkap = preff ** kappa 364 ! WRITE(*,*) 'unskap, cpp, preff :', unskap, cpp, preff 365 DO l = 1, llm 366 DO j=1,jjp1 367 DO i =1, iip1 368 pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap 369 ENDDO 370 ENDDO 371 ENDDO 372 ! 373 ! WRITE(*,*) 'PLS :', pls(10,20,:) 374 ! 375 varname = 'surfgeo' 376 phis(:,:) = 0.0 377 CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 , 378 , jjm ,rlonu,rlatv, interbar ) 379 ! 380 varname = 'u' 381 uvent(:,:,:) = 0.0 382 CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls, 383 . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar ) 384 ! 385 varname = 'v' 386 vvent(:,:,:) = 0.0 387 CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls, 388 . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar ) 389 ! 390 varname = 't' 391 t3d(:,:,:) = 0.0 392 CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, 393 . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar ) 394 ! 395 WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)), 396 . maxval(t3d(:,:,:)) 397 varname = 'tpot' 398 tpot(:,:,:) = 0.0 399 CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, 400 . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar ) 401 ! 402 WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)), 403 . maxval(t3d(:,:,:)) 404 WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)), 405 . maxval(pls(:,:,:)) 406 407 c Calcul de l'humidite a saturation 408 print*,'avant q_sat' 409 call q_sat(llm*jjp1*iip1,t3d,pls,qsat) 410 print*,'apres q_sat' 411 412 WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)), 413 . maxval(qsat(:,:,:)) 414 ! 415 CC WRITE(*,*) 'QSAT :', qsat(10,20,:) 416 ! 417 varname = 'q' 418 qd(:,:,:) = 0.0 419 q3d(:,:,:,:) = 0.0 420 WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)), 421 . maxval(qsat(:,:,:)) 422 CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls, 423 . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar ) 424 q3d(:,:,:,1) = qd(:,:,:) 425 ! 426 427 ! Ozone climatology: 428 if (read_climoz >= 1) call regr_lat_time_climoz(read_climoz) 429 430 varname = 'tsol' 431 ! This line needs to be replaced by a call to restget to get the values in the restart file 432 tsol(:) = 0.0 433 CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0, 434 . jjm, rlonu, rlatv , interbar ) 435 ! 436 WRITE(*,*) 'TSOL construit :' 437 ! WRITE(*,'(48I3)') INT(TSOL(2:klon)-273) 438 ! 439 varname = 'qsol' 440 qsol(:) = 0.0 441 CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0, 442 . jjm, rlonu, rlatv , interbar ) 443 ! 444 varname = 'snow' 445 sn(:) = 0.0 446 CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0, 447 . jjm, rlonu, rlatv , interbar ) 448 ! 449 varname = 'rads' 450 radsol(:) = 0.0 451 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0, 452 . jjm, rlonu, rlatv , interbar ) 453 ! 454 varname = 'rugmer' 455 rugmer(:) = 0.0 456 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0, 457 . jjm, rlonu, rlatv , interbar ) 458 ! 459 ! varname = 'agesno' 460 ! agesno(:) = 0.0 461 ! CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0, 462 ! . jjm, rlonu, rlatv , interbar ) 463 464 varname = 'zmea' 465 zmea(:) = 0.0 466 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0, 467 . jjm, rlonu, rlatv , interbar ) 468 469 varname = 'zstd' 470 zstd(:) = 0.0 471 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0, 472 . jjm, rlonu, rlatv , interbar ) 473 varname = 'zsig' 474 zsig(:) = 0.0 475 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0, 476 . jjm, rlonu, rlatv , interbar ) 477 varname = 'zgam' 478 zgam(:) = 0.0 479 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0, 480 . jjm, rlonu, rlatv , interbar ) 481 varname = 'zthe' 482 zthe(:) = 0.0 483 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0, 484 . jjm, rlonu, rlatv , interbar ) 485 varname = 'zpic' 486 zpic(:) = 0.0 487 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0, 488 . jjm, rlonu, rlatv , interbar ) 489 varname = 'zval' 490 zval(:) = 0.0 491 CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0, 492 . jjm, rlonu, rlatv , interbar ) 493 c 494 cc rugsrel(:) = 0.0 495 cc IF(ok_orodr) THEN 496 cc DO i = 1, iip1* jjp1 497 cc rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. ) 498 cc ENDDO 499 cc ENDIF 500 501 502 C 503 C lecture du fichier glace de terre pour fixer la fraction de terre 504 C et de glace de terre 505 C 506 CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp 507 $ , fid) 508 ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret) 509 ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret) 510 ALLOCATE(dlon_lic(iml_lic), stat=iret) 511 ALLOCATE(dlat_lic(jml_lic), stat=iret) 512 ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret) 513 CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp 514 $ , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) 515 CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp 516 $ , 1, 1, fraclic) 517 CALL flinclo(fid) 518 C 519 C interpolation sur la grille T du modele 520 C 521 WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ', 522 $ iml_lic, jml_lic 523 c 524 C sil les coordonnees sont en degres, on les transforme 525 C 526 IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) ) THEN 527 lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180. 528 ENDIF 529 IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN 530 lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180. 531 ENDIF 532 533 dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1) 534 dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic) 535 C 536 CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic 537 $ ,iim, jjp1, 538 $ rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1)) 539 cx$$$ flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1) 540 flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1) 541 C 542 C passage sur la grille physique 543 C 544 CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, 545 $ pctsrf(1:klon, is_lic)) 546 C adequation avec le maque terre/mer 547 c zmasq(157) = 0. 548 WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA ) 549 pctsrf(1 : klon, is_lic) = 0. 550 END WHERE 551 WHERE (zmasq( 1 : klon) .LT. EPSFRA) 552 pctsrf(1 : klon, is_lic) = 0. 553 END WHERE 554 pctsrf(1 : klon, is_ter) = zmasq(1 : klon) 555 DO ji = 1, klon 556 IF (zmasq(ji) .GT. EPSFRA) THEN 557 IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN 558 pctsrf(ji, is_lic) = zmasq(ji) 559 pctsrf(ji, is_ter) = 0. 560 ELSE 561 pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic) 562 IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN 563 pctsrf(ji,is_ter) = 0. 564 pctsrf(ji, is_lic) = zmasq(ji) 565 ENDIF 566 ENDIF 567 ENDIF 568 END DO 569 C 570 C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0) 571 C 572 pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon)) 573 574 575 WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA) 576 pctsrf(1 : klon, is_oce) = 0. 577 END WHERE 578 579 if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon) 580 581 isst = 0 582 where (pctsrf(2:klon-1,is_oce) >0.) isst = 1 583 C 584 C verif que somme des sous surface = 1 585 C 586 ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf),dim=2))-1.0) 587 $ .GT. EPSFRA) 588 IF (ji .NE. 0) THEN 589 WRITE(*,*) 'pb repartition sous maille pour ',ji,' points' 590 ENDIF 591 592 ! where (pctsrf(1:klon, is_ter) >= .5) 593 ! pctsrf(1:klon, is_ter) = 1. 594 ! pctsrf(1:klon, is_oce) = 0. 595 ! pctsrf(1:klon, is_sic) = 0. 596 ! pctsrf(1:klon, is_lic) = 0. 597 ! zmasq = 1. 598 ! endwhere 599 ! where (pctsrf(1:klon, is_lic) >= .5) 600 ! pctsrf(1:klon, is_ter) = 0. 601 ! pctsrf(1:klon, is_oce) = 0. 602 ! pctsrf(1:klon, is_sic) = 0. 603 ! pctsrf(1:klon, is_lic) = 1. 604 ! zmasq = 1. 605 ! endwhere 606 ! where (pctsrf(1:klon, is_oce) >= .5) 607 ! pctsrf(1:klon, is_ter) = 0. 608 ! pctsrf(1:klon, is_oce) = 1. 609 ! pctsrf(1:klon, is_sic) = 0. 610 ! pctsrf(1:klon, is_lic) = 0. 611 ! zmasq = 0. 612 ! endwhere 613 ! where (pctsrf(1:klon, is_sic) >= .5) 614 ! pctsrf(1:klon, is_ter) = 0. 615 ! pctsrf(1:klon, is_oce) = 0. 616 ! pctsrf(1:klon, is_sic) = 1. 617 ! pctsrf(1:klon, is_lic) = 0. 618 ! zmasq = 0. 619 ! endwhere 620 ! call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) 621 C 622 C verif que somme des sous surface = 1 623 C 624 ! ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 ) 625 ! $ .GT. EPSFRA) 626 ! IF (ji .NE. 0) THEN 627 ! WRITE(*,*) 'pb repartition sous maille pour ',ji,' points' 628 ! ENDIF 629 630 CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d) 631 write(*,*)'zmasq = ' 632 write(*,'(96i1)')nint(zx_tmp_2d) 633 call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) 634 WRITE(*,*) 'MASQUE construit : Masque' 635 WRITE(*,'(97I1)') nINT(masque(:,:)) 636 637 638 639 C Calcul intermediaire 640 c 641 CALL massdair( p3d, masse ) 642 c 643 644 print *,' ALPHAX ',alphax 645 646 DO l = 1, llm 647 DO i = 1, iim 648 xppn(i) = aire( i, 1 ) * masse( i , 1 , l ) 649 xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l ) 650 ENDDO 651 xpn = SUM(xppn)/apoln 652 xps = SUM(xpps)/apols 653 DO i = 1, iip1 654 masse( i , 1 , l ) = xpn 655 masse( i , jjp1 , l ) = xps 656 ENDDO 657 ENDDO 658 q3d(iip1,:,:,:) = q3d(1,:,:,:) 659 phis(iip1,:) = phis(1,:) 660 661 C Ecriture 662 CALL inidissip( lstardis, nitergdiv, nitergrot, niterh , 663 * tetagdiv, tetagrot , tetatemp ) 664 print*,'sortie inidissip' 665 itau = 0 666 itau_dyn = 0 667 itau_phy = 0 668 iday = dayref +itau/day_step 669 time = real(itau-(iday-dayref)*day_step)/day_step 670 c 671 IF(time.GT.1) THEN 672 time = time - 1 673 iday = iday + 1 674 ENDIF 675 day_ref = dayref 676 annee_ref = anneeref 677 678 CALL geopot ( ip1jmp1, tpot , pk , pks, phis , phi ) 679 print*,'sortie geopot' 680 681 CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis , 682 * phi,w, pbaru,pbarv,time+iday-dayref ) 683 print*,'sortie caldyn0' 684 CALL dynredem0("start.nc",dayref,phis) 685 print*,'sortie dynredem0' 686 CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,masse , 687 . psol) 688 print*,'sortie dynredem1' 689 C 690 C Ecriture etat initial physique 691 C 692 write(*,*)'phystep ',dtvr,iphysiq,nbapp_rad 693 phystep = dtvr * FLOAT(iphysiq) 694 radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) ) 695 write(*,*)'phystep =', phystep, radpas 696 cIM : lecture de co2_ppm & solaire ds physiq.def 697 c co2_ppm = 348.0 698 c solaire = 1365.0 699 700 c 701 c Initialisation 702 c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs 703 c 704 ftsol(:,is_ter) = tsol 705 ftsol(:,is_lic) = tsol 706 ftsol(:,is_oce) = tsol 707 ftsol(:,is_sic) = tsol 708 snsrf(:,is_ter) = sn 709 snsrf(:,is_lic) = sn 710 snsrf(:,is_oce) = sn 711 snsrf(:,is_sic) = sn 712 falb1(:,is_ter) = 0.08 713 falb1(:,is_lic) = 0.6 714 falb1(:,is_oce) = 0.5 715 falb1(:,is_sic) = 0.6 716 falb2 = falb1 717 evap(:,:) = 0. 718 qsolsrf(:,is_ter) = 150 719 qsolsrf(:,is_lic) = 150 720 qsolsrf(:,is_oce) = 150. 721 qsolsrf(:,is_sic) = 150. 722 do i = 1, nbsrf 723 do j = 1, nsoilmx 724 tsoil(:,j,i) = tsol 725 enddo 726 enddo 727 rain_fall = 0.; snow_fall = 0. 728 solsw = 165. 729 sollw = -53. 730 t_ancien = 273.15 731 q_ancien = 0. 732 agesno = 0. 733 c 734 frugs(1:klon,is_oce) = rugmer(1:klon) 735 frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0) 736 frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0) 737 frugs(1:klon,is_sic) = 0.001 738 fder = 0.0 739 clwcon = 0.0 740 rnebcon = 0.0 741 ratqs = 0.0 742 run_off_lic_0 = 0.0 743 rugoro = 0.0 744 745 c 746 c Avant l'appel a phyredem, on initialize les modules de surface 747 c avec les valeurs qui vont etre ecrit dans startphy.nc 748 c 749 dummy = 1.0 750 pbl_tke(:,:,:) = 1.e-8 751 zmax0(:) = 40. 752 f0(:) = 1.e-5 753 ema_work1(:,:) = 0. 754 ema_work2(:,:) = 0. 755 wake_deltat(:,:) = 0. 756 wake_deltaq(:,:) = 0. 757 wake_s(:) = 0. 758 wake_cstar(:) = 0. 759 wake_fip(:) = 0. 760 761 call fonte_neige_init(run_off_lic_0) 762 call pbl_surface_init(qsol, fder, snsrf, qsolsrf, 763 $ evap, frugs, agesno, tsoil) 764 765 call phyredem("startphy.nc") 766 767 768 769 C Sortie Visu pour les champs dynamiques 770 cc if (1.eq.0 ) then 771 cc print*,'sortie visu' 772 cc time_step = 1. 773 cc t_ops = 2. 774 cc t_wrt = 2. 775 cc itau = 2. 776 cc visu_file='Etat0_visu.nc' 777 cc CALL initdynav(visu_file,dayref,anneeref,time_step, 778 cc . t_ops, t_wrt, visuid) 779 cc CALL writedynav(visuid, itau,vvent , 780 cc . uvent,tpot,pk,phi,q3d,masse,psol,phis) 781 cc else 782 print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0' 783 cc endif 784 print*,'entree histclo' 785 CALL histclo 137 co2_ppm0 = co2_ppm 138 139 dtvr = daysec/FLOAT(day_step) 140 WRITE(lunout,*)'dtvr',dtvr 141 142 CALL iniconst() 143 CALL inigeom() 144 145 !--- Initializations for tracers 146 CALL infotrac_init 147 ALLOCATE(q3d(iip1,jjp1,llm,nqtot)) 148 149 CALL inifilr() 150 CALL phys_state_var_init(read_climoz) 151 152 rlat(1) = ASIN(1.0) 153 DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j); END DO 154 rlat(klon) = - ASIN(1.0) 155 rlat(:)=rlat(:)*(180.0/pi) 156 157 rlon(1) = 0.0 158 DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO 159 rlon(klon) = 0.0 160 rlon(:)=rlon(:)*(180.0/pi) 161 162 ! For a coupled simulation, the ocean mask from ocean model is used to compute 163 ! the weights an to insure ocean fractions are the same for atmosphere and ocean 164 ! Otherwise, mask is created using Relief file. 165 166 WRITE(lunout,*)'Essai de lecture masque ocean' 167 iret = NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a) 168 IF(iret/=NF90_NOERR) THEN 169 WRITE(lunout,*)'ATTENTION!! pas de fichier o2a.nc trouve' 170 WRITE(lunout,*)'Run force' 171 x='masque' 172 masque(:,:)=0.0 173 CALL startget(x, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm, rlonu, & 174 rlatv, ib) 175 WRITE(lunout,*)'MASQUE construit : Masque' 176 WRITE(lunout,'(97I1)') nINT(masque) 177 CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq) 178 WHERE( zmasq(:)<EPSFRA) zmasq(:)=0. 179 WHERE(1.-zmasq(:)<EPSFRA) zmasq(:)=1. 180 ELSE 181 couple=.true. 182 iret=NF90_CLOSE(nid_o2a) 183 CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a) 184 IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN 185 WRITE(lunout,*)'Dimensions non compatibles pour masque ocean' 186 WRITE(lunout,*)'iim = ',iim,' iml_omask = ',iml_omask 187 WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask 188 CALL abort_gcm('etat0_netcdf','Dimensions non compatibles pour masque oc& 189 &ean',1) 190 END IF 191 ALLOCATE( ocemask(iml_omask,jml_omask), ocetmp(iml_omask,jml_omask)) 192 ALLOCATE(lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask)) 193 ALLOCATE(dlon_omask(iml_omask),dlat_omask(jml_omask)) 194 CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, & 195 lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid) 196 CALL flinget(fid,'OceMask',iml_omask,jml_omask,llm_tmp,ttm_tmp,1,1,ocetmp) 197 CALL flinclo(fid) 198 dlon_omask(:)=lon_omask(:,1) 199 dlat_omask(:)=lat_omask(1,:) 200 ocemask=ocetmp 201 IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN 202 DO j=1,jml_omask 203 ocemask(:,j) = ocetmp(:,jml_omask-j+1) 204 END DO 205 END IF 206 ALLOCATE( ocemask(iml_omask,jml_omask), ocetmp(iml_omask,jml_omask)) 207 ALLOCATE( lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask)) 208 ALLOCATE(dlon_omask(iml_omask), dlat_omask(jml_omask)) 209 CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, lon_omask, & 210 lat_omask, lev, ttm_tmp, itaul, date, dt, fid) 211 CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, ttm_tmp, & 212 1, 1, ocetmp) 213 CALL flinclo(fid) 214 dlon_omask(1:iml_omask) = lon_omask(1:iml_omask,1) 215 dlat_omask(1:jml_omask) = lat_omask(1,1:jml_omask) 216 ocemask = ocetmp 217 IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN 218 DO j=1,jml_omask 219 ocemask(:,j) = ocetmp(:,jml_omask-j+1) 220 END DO 221 END IF 222 ! 223 ! Ocean mask to physical grid 224 !******************************************************************************* 225 WRITE(lunout,*)'ocemask ' 226 WRITE(fmt,"(i4,'i1)')")iml_omask ; fmt='('//ADJUSTL(fmt) 227 WRITE(lunout,fmt)int(ocemask) 228 ocemask_fi(1)=ocemask(1,1) 229 DO j=2,jjm; ocemask_fi((j-2)*iim+1:iim+1)=ocemask(1:iim,j); END DO 230 ocemask_fi(klon)=ocemask(1,jjp1) 231 zmasq=1.-ocemask_fi 232 END IF 233 234 CALL gr_fi_dyn(1,klon,iip1,jjp1,zmasq,masque) 235 236 ! The startget calls need to be replaced by a call to restget to get the 237 ! values in the restart file 238 x = 'relief'; orog(:,:) = 0.0 239 CALL startget(x,iip1,jjp1,rlonv,rlatu, orog, 0.0,jjm,rlonu,rlatv,ib,& 240 masque) 241 242 x = 'rugosite'; rugo(:,:) = 0.0 243 CALL startget(x,iip1,jjp1,rlonv,rlatu, rugo, 0.0,jjm, rlonu,rlatv,ib) 244 ! WRITE(lunout,'(49I1)') INT(orog(:,:)*10) 245 ! WRITE(lunout,'(49I1)') INT(rugo(:,:)*10) 246 247 ! Sub-surfaces initialization 248 !******************************************************************************* 249 pctsrf=0. 250 x = 'psol'; psol(:,:) = 0.0 251 CALL startget(x,iip1,jjp1,rlonv,rlatu,psol,0.0,jjm,rlonu,rlatv,ib) 252 ! WRITE(lunout,*) 'PSOL :', psol(10,20) 253 ! WRITE(lunout,*) ap(:), bp(:) 254 255 ! Mid-levels pressure computation 256 !******************************************************************************* 257 CALL pression(ip1jmp1, ap, bp, psol, p3d) 258 CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y) 259 pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa) 260 ! WRITE(lunout,*) 'P3D :', p3d(10,20,:) 261 ! WRITE(lunout,*) 'PK:', pk(10,20,:) 262 ! WRITE(lunout,*) 'PLS :', pls(10,20,:) 263 264 x = 'surfgeo'; phis(:,:) = 0.0 265 CALL startget(x,iip1,jjp1,rlonv,rlatu, phis, 0.0,jjm, rlonu,rlatv,ib) 266 267 x = 'u'; uvent(:,:,:) = 0.0 268 CALL startget(x,iip1,jjp1,rlonu,rlatu,llm,pls,y,uvent,0.0,jjm, rlonv,rlatv,ib) 269 270 x = 'v'; vvent(:,:,:) = 0.0 271 CALL startget(x,iip1,jjm, rlonv,rlatv,llm,pls,y,vvent,0.0,jjp1,rlonu,rlatu,ib) 272 273 x = 't'; t3d(:,:,:) = 0.0 274 CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,y,t3d, 0.0,jjm, rlonu,rlatv,ib) 275 276 x = 'tpot'; tpot(:,:,:) = 0.0 277 CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,pk,tpot,0.0,jjm, rlonu,rlatv,ib) 278 279 WRITE(lunout,*) 'T3D min,max:',minval(t3d(:,:,:)),maxval(t3d(:,:,:)) 280 WRITE(lunout,*) 'PLS min,max:',minval(pls(:,:,:)),maxval(pls(:,:,:)) 281 282 ! Humidity at saturation computation 283 !******************************************************************************* 284 WRITE(lunout,*) 'avant q_sat' 285 CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat) 286 WRITE(lunout,*) 'apres q_sat' 287 WRITE(lunout,*) 'QSAT min,max:',minval(qsat(:,:,:)),maxval(qsat(:,:,:)) 288 ! WRITE(lunout,*) 'QSAT :',qsat(10,20,:) 289 290 x = 'q'; qd (:,:,:) = 0.0 291 CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,qsat,qd,0.0,jjm, rlonu,rlatv,ib) 292 q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:) 293 294 !--- OZONE CLIMATOLOGY 295 IF(read_climoz>=1) CALL regr_lat_time_climoz(read_climoz) 296 297 x = 'tsol'; tsol(:) = 0.0 298 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,tsol, 0.0,jjm,rlonu,rlatv,ib) 299 300 x = 'qsol'; qsol(:) = 0.0 301 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,qsol, 0.0,jjm,rlonu,rlatv,ib) 302 303 x = 'snow'; sn(:) = 0.0 304 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,sn, 0.0,jjm,rlonu,rlatv,ib) 305 306 x = 'rads'; radsol(:) = 0.0 307 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,jjm,rlonu,rlatv,ib) 308 309 x = 'rugmer'; rugmer(:) = 0.0 310 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,jjm,rlonu,rlatv,ib) 311 312 x = 'zmea'; zmea(:) = 0.0 313 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zmea, 0.0,jjm,rlonu,rlatv,ib) 314 315 x = 'zstd'; zstd(:) = 0.0 316 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zstd, 0.0,jjm,rlonu,rlatv,ib) 317 318 x = 'zsig'; zsig(:) = 0.0 319 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zsig, 0.0,jjm,rlonu,rlatv,ib) 320 321 x = 'zgam'; zgam(:) = 0.0 322 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zgam, 0.0,jjm,rlonu,rlatv,ib) 323 324 x = 'zthe'; zthe(:) = 0.0 325 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zthe, 0.0,jjm,rlonu,rlatv,ib) 326 327 x = 'zpic'; zpic(:) = 0.0 328 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zpic, 0.0,jjm,rlonu,rlatv,ib) 329 330 x = 'zval'; zval(:) = 0.0 331 CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zval, 0.0,jjm,rlonu,rlatv,ib) 332 333 ! WRITE(lunout,'(48I3)') 'TSOL :', INT(tsol(2:klon)-273) 334 335 ! Soil ice file reading for soil fraction and soil ice fraction 336 !******************************************************************************* 337 CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid) 338 ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic)) 339 ALLOCATE(dlat_lic(jml_lic), dlon_lic(iml_lic)) 340 ALLOCATE( fraclic(iml_lic,jml_lic)) 341 CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp, & 342 lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) 343 CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic) 344 CALL flinclo(fid) 345 346 ! Interpolation on model T-grid 347 !******************************************************************************* 348 WRITE(lunout,*)'dimensions de landice iml_lic, jml_lic : ',iml_lic,jml_lic 349 ! conversion if coordinates are in degrees 350 IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. 351 IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180. 352 dlon_lic(:)=lon_lic(:,1) 353 dlat_lic(:)=lat_lic(1,:) 354 CALL grille_m( iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic, iim,jjp1, & 355 rlonv, rlatu, flic_tmp(1:iim,:) ) 356 flic_tmp(iip1,:)=flic_tmp(1,:) 357 358 !--- To the physical grid 359 CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic)) 360 361 !--- Adequation with soil/sea mask 362 WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 363 WHERE(zmasq(:)<EPSFRA) pctsrf(:,is_lic)=0. 364 pctsrf(:,is_ter)=zmasq(:) 365 DO ji=1,klon 366 IF(zmasq(ji)>EPSFRA) THEN 367 IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN 368 pctsrf(ji,is_lic)=zmasq(ji) 369 pctsrf(ji,is_ter)=0. 370 ELSE 371 pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic) 372 IF(pctsrf(ji,is_ter)<EPSFRA) THEN 373 pctsrf(ji,is_ter)=0. 374 pctsrf(ji,is_lic)=zmasq(ji) 375 END IF 376 END IF 377 END IF 378 END DO 379 380 ! sub-surface ocean and sea ice (sea ice set to zero for start) 381 !******************************************************************************* 382 pctsrf(:,is_oce)=(1.-zmasq(:)) 383 WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0. 384 IF(couple) pctsrf(:,is_oce)=ocemask_fi(:) 385 isst=0 386 WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1 387 388 ! It is checked that the sub-surfaces sum is equal to 1 389 !******************************************************************************* 390 ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA) 391 IF(ji/=0) WRITE(lunout,*) 'pb repartition sous maille pour ',ji,' points' 392 CALL gr_fi_ecrit(1, klon, iim, jjp1, zmasq, zx_tmp_2d) 393 ! WRITE(fmt,"(i3,')')")iim; fmt='(i'//ADJUSTL(fmt) 394 ! WRITE(lunout,*)'zmasq = ' 395 ! WRITE(lunout,TRIM(fmt))NINT(zx_tmp_2d) 396 CALL gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque) 397 WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt) 398 WRITE(lunout,*) 'MASQUE construit : Masque' 399 WRITE(lunout,TRIM(fmt))NINT(masque(:,:)) 400 401 ! Intermediate computation 402 !******************************************************************************* 403 CALL massdair(p3d,masse) 404 WRITE(lunout,*)' ALPHAX ',alphax 405 DO l=1,llm 406 xppn(:)=aire(1:iim,1 )*masse(1:iim,1 ,l) 407 xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l) 408 xpn=SUM(xppn)/apoln 409 xps=SUM(xpps)/apols 410 masse(:,1 ,l)=xpn 411 masse(:,jjp1,l)=xps 412 END DO 413 q3d(iip1,:,:,:)=q3d(1,:,:,:) 414 phis(iip1,:) = phis(1,:) 415 416 IF(.NOT.letat0) RETURN 417 418 ! Writing 419 !******************************************************************************* 420 CALL inidissip(lstardis,nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,tetatemp) 421 WRITE(lunout,*)'sortie inidissip' 422 itau=0 423 itau_dyn=0 424 itau_phy=0 425 iday=dayref+itau/day_step 426 time=FLOAT(itau-(iday-dayref)*day_step)/day_step 427 IF(time>1.) THEN 428 time=time-1 429 iday=iday+1 430 END IF 431 day_ref=dayref 432 annee_ref=anneeref 433 434 CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi ) 435 WRITE(lunout,*)'sortie geopot' 436 437 CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, & 438 phi, w, pbaru, pbarv, time+iday-dayref) 439 WRITE(lunout,*)'sortie caldyn0' 440 CALL dynredem0( "start.nc", dayref, phis) 441 WRITE(lunout,*)'sortie dynredem0' 442 CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) 443 WRITE(lunout,*)'sortie dynredem1' 444 445 ! Physical initial state writting 446 !******************************************************************************* 447 WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad 448 phystep = dtvr * FLOAT(iphysiq) 449 radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) ) 450 WRITE(lunout,*)'phystep =', phystep, radpas 451 452 ! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs 453 !******************************************************************************* 454 DO i=1,nbsrf; ftsol(:,i) = tsol; END DO 455 DO i=1,nbsrf; snsrf(:,i) = sn; END DO 456 falb1(:,is_ter) = 0.08; falb1(:,is_lic) = 0.6 457 falb1(:,is_oce) = 0.5; falb1(:,is_sic) = 0.6 458 falb2 = falb1 459 evap(:,:) = 0. 460 DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO 461 DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO 462 rain_fall = 0.; snow_fall = 0. 463 solsw = 165.; sollw = -53. 464 t_ancien = 273.15 465 q_ancien = 0. 466 agesno = 0. 467 frugs(:,is_oce) = rugmer(:) 468 frugs(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) 469 frugs(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) 470 frugs(:,is_sic) = 0.001 471 fder = 0.0 472 clwcon = 0.0 473 rnebcon = 0.0 474 ratqs = 0.0 475 run_off_lic_0 = 0.0 476 rugoro = 0.0 477 478 ! Before phyredem calling, surface modules and values to be saved in startphy.nc 479 ! are initialized 480 !******************************************************************************* 481 dummy = 1.0 482 pbl_tke(:,:,:) = 1.e-8 483 zmax0(:) = 40. 484 f0(:) = 1.e-5 485 ema_work1(:,:) = 0. 486 ema_work2(:,:) = 0. 487 wake_deltat(:,:) = 0. 488 wake_deltaq(:,:) = 0. 489 wake_s(:) = 0. 490 wake_cstar(:) = 0. 491 wake_fip(:) = 0. 492 CALL fonte_neige_init(run_off_lic_0) 493 CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil ) 494 CALL phyredem( "startphy.nc" ) 495 496 WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0' 497 WRITE(lunout,*)'entree histclo' 498 CALL histclo() 786 499 787 500 #endif 788 501 !#endif of #ifdef CPP_EARTH 789 RETURN 790 ! 791 END SUBROUTINE etat0_netcdf 502 RETURN 503 504 END SUBROUTINE etat0_netcdf 505 ! 506 !-------------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
