--- trunk/libf/dyn3d/etat0.f90 2008/03/03 16:32:04 5 +++ trunk/libf/dyn3d/etat0.f90 2008/08/07 12:29:13 18 @@ -6,6 +6,8 @@ IMPLICIT NONE REAL pctsrf(klon, nbsrf) + ! ("pctsrf(i, :)" is the composition of the surface at horizontal + ! position "i") private nbsrf, klon @@ -15,11 +17,11 @@ ! From "etat0_netcdf.F", version 1.3 2005/05/25 13:10:09 - ! This subroutine creates "masque". + ! This subroutine creates "mask". USE ioipsl, only: flinget, flinclo, flinopen_nozoom, flininfo, histclo - USE start_init_orog_m, only: start_init_orog, masque, phis + USE start_init_orog_m, only: start_init_orog, mask, phis use start_init_phys_m, only: qsol_2d use startdyn, only: start_inter_3d, start_init_dyn use dimens_m, only: iim, jjm, llm, nqmx @@ -36,16 +38,16 @@ use serre, only: alphax use dimsoil, only: nsoilmx use temps, only: itau_dyn, itau_phy, annee_ref, day_ref, dt - use clesphys, only: ok_orodr, nbapp_rad use grid_atob, only: grille_m use grid_change, only: init_dyn_phy, dyn_phy use q_sat_m, only: q_sat use exner_hyb_m, only: exner_hyb - use regr_coefoz_m, only: regr_coefoz - use advtrac_m, only: iniadvtrac - use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf90_nowrite, & - nf90_get_var, handle_err - use pressure_m, only: pls, p3d + use iniadvtrac_m, only: iniadvtrac + use pressure_var, only: pls, p3d + use dynredem0_m, only: dynredem0 + use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz + use regr_pr_o3_m, only: regr_pr_o3 + use phyredem_m, only: phyredem ! Variables local to the procedure: @@ -81,7 +83,6 @@ REAL zsig(klon), zgam(klon) REAL zthe(klon) REAL zpic(klon), zval(klon) - REAL rugsrel(klon) REAL t_ancien(klon, llm), q_ancien(klon, llm) ! REAL run_off_lic_0(klon) real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm) @@ -104,8 +105,6 @@ REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) REAL w(ip1jmp1, llm) REAL phystep - INTEGER radpas - integer ncid, varid, ncerr, month !--------------------------------- @@ -132,15 +131,15 @@ lonfi(klon) = 0. call start_init_orog(relief, zstd_2d, zsig_2d, zgam_2d, zthe_2d, zpic_2d, & - zval_2d) ! also compute "masque" and "phis" + zval_2d) ! also compute "mask" and "phis" call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points - zmasq = pack(masque, dyn_phy) + zmasq = pack(mask, dyn_phy) PRINT *, 'Masque construit' CALL start_init_dyn(tsol_2d, psol) ! also compute "qsol_2d" ! Compute pressure on intermediate levels: - forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol(:, :) + forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol CALL exner_hyb(psol, p3d, pks, pk) IF (MINVAL(pk) == MAXVAL(pk)) stop '"pk" should not be constant' @@ -187,46 +186,26 @@ if (nqmx >= 5) then ! Ozone: - - ! Compute ozone parameters on the LMDZ grid: - call regr_coefoz - - ! Find the month containing the day number "dayref": - month = (dayref - 1) / 30 + 1 - print *, "month = ", month - - call nf95_open("coefoz_LMDZ.nc", nf90_nowrite, ncid) - - ! Get data at the right month from the input file: - call nf95_inq_varid(ncid, "r_Mob", varid) - ncerr = nf90_get_var(ncid, varid, q3d(:, :, :, 5), & - start=(/1, 1, 1, month/)) - call handle_err("nf90_get_var r_Mob", ncerr) - - call nf95_close(ncid) - ! Latitudes are in increasing order in the input file while - ! "rlatu" is in decreasing order so we need to invert order. Also, we - ! compute mass fraction from mole fraction: - q3d(:, :, :, 5) = q3d(:, jjm+1:1:-1, :, 5) * 48. / 29. + call regr_lat_time_coefoz + call regr_pr_o3(q3d(:, :, :, 5)) + ! Convert from mole fraction to mass fraction: + q3d(:, :, :, 5) = q3d(:, :, :, 5) * 48. / 29. end if - tsol(:) = pack(tsol_2d, dyn_phy) - qsol(:) = pack(qsol_2d, dyn_phy) - sn(:) = 0. ! snow - radsol(:) = 0. - tslab(:) = 0. ! IM "slab" ocean - seaice(:) = 0. - rugmer(:) = 0.001 - zmea(:) = pack(relief, dyn_phy) - zstd(:) = pack(zstd_2d, dyn_phy) - zsig(:) = pack(zsig_2d, dyn_phy) - zgam(:) = pack(zgam_2d, dyn_phy) - zthe(:) = pack(zthe_2d, dyn_phy) - zpic(:) = pack(zpic_2d, dyn_phy) - zval(:) = pack(zval_2d, dyn_phy) - - rugsrel(:) = 0. - IF (ok_orodr) rugsrel(:) = MAX(1.e-05, zstd(:) * zsig(:) / 2) + tsol = pack(tsol_2d, dyn_phy) + qsol = pack(qsol_2d, dyn_phy) + sn = 0. ! snow + radsol = 0. + tslab = 0. ! IM "slab" ocean + seaice = 0. + rugmer = 0.001 + zmea = pack(relief, dyn_phy) + zstd = pack(zstd_2d, dyn_phy) + zsig = pack(zsig_2d, dyn_phy) + zgam = pack(zgam_2d, dyn_phy) + zthe = pack(zthe_2d, dyn_phy) + zpic = pack(zpic_2d, dyn_phy) + zval = pack(zval_2d, dyn_phy) ! On initialise les sous-surfaces: ! Lecture du fichier glace de terre pour fixer la fraction de terre @@ -251,48 +230,50 @@ print *, "jml_lic = ", jml_lic ! Si les coordonnées sont en degrés, on les transforme : - IF (MAXVAL( lon_lic(:, :) ) > pi) THEN - lon_lic(:, :) = lon_lic(:, :) * pi / 180. + IF (MAXVAL( lon_lic ) > pi) THEN + lon_lic = lon_lic * pi / 180. ENDIF - IF (maxval( lat_lic(:, :) ) > pi) THEN - lat_lic(:, :) = lat_lic(:, :) * pi/ 180. + IF (maxval( lat_lic ) > pi) THEN + lat_lic = lat_lic * pi/ 180. ENDIF - dlon_lic(:) = lon_lic(:, 1) - dlat_lic(:) = lat_lic(1, :) + dlon_lic = lon_lic(:, 1) + dlat_lic = lat_lic(1, :) flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), & rlatu) flic_tmp(iim + 1, :) = flic_tmp(1, :) ! Passage sur la grille physique - pctsrf(:, :)=0. + pctsrf = 0. pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy) ! Adéquation avec le maque terre/mer WHERE (pctsrf(:, is_lic) < EPSFRA ) pctsrf(:, is_lic) = 0. - WHERE (zmasq(:) < EPSFRA) pctsrf(:, is_lic) = 0. - pctsrf(:, is_ter) = zmasq(:) - where (zmasq(:) > EPSFRA) - where (pctsrf(:, is_lic) >= zmasq(:)) - pctsrf(:, is_lic) = zmasq(:) + WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0. + pctsrf(:, is_ter) = zmasq + where (zmasq > EPSFRA) + where (pctsrf(:, is_lic) >= zmasq) + pctsrf(:, is_lic) = zmasq pctsrf(:, is_ter) = 0. elsewhere - pctsrf(:, is_ter) = zmasq(:) - pctsrf(:, is_lic) + pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic) where (pctsrf(:, is_ter) < EPSFRA) pctsrf(:, is_ter) = 0. - pctsrf(:, is_lic) = zmasq(:) + pctsrf(:, is_lic) = zmasq end where end where end where ! Sous-surface océan et glace de mer (pour démarrer on met glace ! de mer à 0) : - pctsrf(:, is_oce) = 1. - zmasq(:) + pctsrf(:, is_oce) = 1. - zmasq WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0. ! Vérification que somme des sous-surfaces vaut 1: - ji = count(abs(sum(pctsrf(:, :), dim = 2) - 1. ) > EPSFRA) - IF (ji /= 0) PRINT *, 'Problème répartition sous maille pour', ji, 'points' + ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA) + IF (ji /= 0) then + PRINT *, 'Problème répartition sous maille pour ', ji, 'points' + end IF ! Calcul intermédiaire: CALL massdair(p3d, masse) @@ -324,11 +305,8 @@ ! Ecriture état initial physique: print *, 'dtvr = ', dtvr print *, "iphysiq = ", iphysiq - print *, "nbapp_rad = ", nbapp_rad phystep = dtvr * REAL(iphysiq) - radpas = NINT (86400./phystep/ nbapp_rad) print *, 'phystep = ', phystep - print *, "radpas = ", radpas ! Initialisations : tsolsrf(:, is_ter) = tsol @@ -344,7 +322,7 @@ albe(:, is_oce) = 0.5 albe(:, is_sic) = 0.6 alblw = albe - evap(:, :) = 0. + evap = 0. qsolsrf(:, is_ter) = 150. qsolsrf(:, is_lic) = 150. qsolsrf(:, is_oce) = 150. @@ -358,12 +336,12 @@ q_ancien = 0. agesno = 0. !IM "slab" ocean - tslab(:) = tsolsrf(:, is_oce) + tslab = tsolsrf(:, is_oce) seaice = 0. - frugs(:, is_oce) = rugmer(:) - frugs(:, is_ter) = MAX(1.e-05, zstd(:) * zsig(:) / 2) - frugs(:, is_lic) = MAX(1.e-05, zstd(:) * zsig(:) / 2) + frugs(:, is_oce) = rugmer + frugs(:, is_ter) = MAX(1.e-05, zstd * zsig / 2) + frugs(:, is_lic) = MAX(1.e-05, zstd * zsig / 2) frugs(:, is_sic) = 0.001 fder = 0. clwcon = 0. @@ -371,10 +349,10 @@ ratqs = 0. run_off_lic_0 = 0. - call phyredem("startphy.nc", phystep, radpas, latfi, lonfi, pctsrf, & + call phyredem("startphy.nc", latfi, lonfi, pctsrf, & tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, & evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, & - agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, & + agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, & t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0) CALL histclo